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Abstract 



In the traveling salesman problem, one must find the length of the shortest closed 
tour visiting given "cities" . We study the stochastic version of the problem, taking the 
locations of cities and the distances separating them to be random variables drawn from an 
ensemble. Wc consider first the ensemble where cities are placed in Euclidean space. We 
investigate how the optimum tour length scales with number of cities and with number 
of spatial dimensions. We then examine the analytical theory behind the random link 
ensemble, where distances between cities are independent random variables. Finally, we 
look at the related geometric issue of nearest neighbor distances, and find some remarkable 
universalities. 
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Chapter I 



Introduction 



Consider A'^ sites ("cities") in a space, and the distances connecting each pair of cities. Now find 
the shortest possible path ("tour") visiting each city exactly once, and returning to its starting 
point. This is the traveling salesman problem (tsp), the most classic and best known problem 
in the field of combinatorial optimization (co). 

The TSP is an extremely simple problem to state; deceptively so, because it is also an 
extremely difficult one to solve. The most naive algorithm for finding the optimum tour would 
have to consider all (N — l)!/2 possible tours. Working this way, the fastest computer on 
the market today would require more time than the current age of the universe to solve a 
case with 27 cities. Of course, far more sophisticated algorithms exist, but even the best of 
these is sharply limited by the inherent complexity of the TSP. Computer scientists classify 
the problem as being NP-hard: roughly speaking this means that, most likely, there exists no 
algorithm that can consistently find the optimal tour in an amount of time polynomial in N. The 
statement provides a technical standard for "hardness" , saying that if a polynomial algorithm 
could ever be found for the tsp, one would exist for a whole class of other co problems as well 
(Garey & Johnson, 1979, see also discussion in Appendix A). For this reason, and because 
of how straightforward the TSP is to state, it has served as a prototype for the study of difficult 
combinatorial optimization problems in general. 

The challenge of making progress on the TSP despite its complexity has carried its appeal 
well beyond its initial boundaries. Research on the TSP and related CO issues has become 
a truly interdisciplinary effort, uniting computer scientists, pure mathematicians, operations 
researchers, and more recently, theoretical physicists.^ In physics, the TSP has been of particular 
interest in its stochastic form, where the configuration ( "instance" ) of the N cities is randomly 
chosen from an ensemble (see Figure I-l for an example of an iV = 24 instance and its optimal 
tour). Starting in the 1980s, Kirkpatrick, Mczard, Parisi and other theoreticians began applying 
tools from statistical mechanics to the stochastic tsp, noting broad similarities between CO 
problems and the emerging field of disordered systems. The hope was, first of all, that numerical 
methods developed in the context of disordered systems could be adapted to the TSP, and 
second of all, that the analytical formalism developed for these other problems could provide a 
theoretical framework for understanding the TSP — and ultimately, could provide exact results 
as well. 

The generic example of a disordered system is the spin glass. A spin glass may be thought of 
as a ferromagnetic Ising model — the fundamental lattice model that statistical physicists use 

^Biology, in the guise of genetic algorithms and DNA computing, is but the latest field to arrive on the scene. 



7 



Chap. I: Introduction 



Figure I— 1: An elementary instance of tlie traveling salesman problem, and its optimal tour. 



to describe magnetism — but modified so that the interaction strengths between neighboring 
lattice sites are random quantities. If, as in the case of the spin glass, they are allowed to be 
negative as well as positive (antiferromagnetic as well as ferromagnetic), we have a frustrated 
disordered system. Take a fixed instance of these random variables ("quenched disorder"). 
The total energy of the system, involving a sum over all pairs of neighboring sites, is then a 
complicated function of the particular state of the system (meaning the binary values of the 
spins located at each of N sites), and may be represented schematically by a "landscape" of 
hills and valleys (see Figure 1-2). We shall concern ourselves with what statistical physicists 
think of as the zero-temperature case, where we must find the lowest-energy state, and hence the 
global minimum in this landscape. If we are to avoid actually enumerating all possible states 
— a task which increases exponentially with the number of sites — the problem of finding the 
global minimum becomes highly nontrivial since there are numerous local minima in which a 
minimizing search procedure could get trapped. 

Let us now examine the TSP using this picture. The system's "energy" is the tour length, 
and is to be minimized. The quenched "interaction strengths" are the distances between cities. 
The energy is once again a complicated function of the state of the system, where in this case 
"state" refers to the particular tour chosen. Optimizing the tour is precisely equivalent to finding 
the system's ground state, and thus the global minimum in Figure 1-2. If wc look at it in this 
way, the same approximation methods used to find spin glass ground states may be used to 
find optimal TSP tours; simulated annealing, developed independently by KiRKPATRiCK (1984) 
and Cerny (1985), is the most famous example. At the same time, if we look not at a single 
instance but rather at the stochastic ensemble of instances, analytical approaches giving exact 
results for spin glass problems can be adapted to provide predictions for the TSP that avoid 
algorithmic procedures altogether. 
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The foregoing analogy can be broadened: disordered systems and CO problems display many 
similar properties, often making methods developed in one field applicable to the other. A 
common feature in combinatorial optimization is the existence of a large number of quantities 
interacting in a complicated manner, where each quantity exerts an influence on the others — 
and these influences are often contradictory. Describing such systems on a microscopic level is 
generally so difficult as to be impossible. A method better suited to these problems is that of 
statistical physics, where one contents oneself with studying the broader picture, reducing the 
number of parameters considerably by resorting to a probabilistic approach. One thus obtains 
the macroscopic quantities of interest, to reasonable accuracy, with far less effort. 

Let us express these concepts more explicitly. A CO problem such as the TSP may be 
formulated as follows. We have an objective function representing "cost" — the tour length, for 
the TSP — that must be minimized. This cost may be expressed as 

Liiriij}) = ^kjnij, (I.l) 

where i and j go from 1 up to AT (the total system "size"), the {lij} are the parameters of the 
problem, and the {riij} are the variables whose values must be adjusted so as to minimize L. A 
certain set of constraints determines the allowed values for the {riij}. In the case of the TSP, for 
example, L is simply the tour length to be minimized, {hj} is the matrix of distances between 
city i and city j, and {riij} is a binary matrix expressing whether the link between i and j is 
used on the tour. The constraints on the variables are thus that: (i) each link is either used once 
or not at all (Vi, j, riij £ {0, 1} and Vi, na = 0); (ii) each city has one link going into it and one 
link coming out from it (Vj, X^jn^ = 1 and Vi, ^ ■ riij = 1); and (iii) the {riij} are such that 



Energy 




Global minimum 

States »■ 

Figure 1—2: Schematic energy "landscape", over state space, of a frustrated disordered system. 
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the links form exactly one loop (Vii, 3i2, ...,iN ■ rii^i^ = rii^i^ = ■■■ = nij^_^ij^ = riij^i^ = 1).^ 

To put this into the language of statistical physics, we think of L as the system's energy; 
each contribution lijUij is then the interaction energy associated with the pair of sites i and j. 
We may write down the partition function: 

Z=Y^ e-^«"*^»/^, (1.2) 

{Uij} 

where the sum is taken over the ensemble of matrices {riij } that satisfy the constraints we have 
posed, and T is the physical "temperature" of the system. In the limit T ^ the free energy 
F = —TlnZ will be dominated by the ground-state contribution, that is, the particular state 
(value of the matrix {riij}) that minimizes L. When we speak of the T = solution, we are thus 
speaking of the global optimum. 

Rephrasing the problem in this manner does not appear to help us find the exact state 
minimizing L for a given set of parameters {hj} — after all, when A'^ is large and we are dealing 
with matrices of size N x N, there is still an enormous number of parameters. However, we may 
instead consider the {kj} as being chosen randomly from an ensemble, and ask for properties 
of the distribution of L, such as its mean over the ensemble. Analytical tools from statistical 
physics, and particularly the approximations that are typically used, then become relevant and 
useful. The stochastic TSP is a prime candidate for this sort of treatment. 

Two main versions of the stochastic tsp have been discussed in the literature, and form the 
subject of the majority of this thesis. The first and more classic version is the Euclidean TSP (also 
known as the random point TSP), where the N cities are placed randomly and independently 
in d-dimensional Euclidean space, and the distances between cities arc given by the Euclidean 
metric. The second is the random link TSP, largely developed within the context of disordered 
systems. In the random link TSP, it is not the positions of the cities but rather the lengths 
separating them that are independent and identically distributed random variables. We speak 
of "lengths" because there is no physical "distance" here, nor geometry of any sort, apart from 
a symmetry requiring that the length be equal to Iji between any two cities i and j.^ The 
great advantage of the random link tsp over the Euclidean (random point) tsp is that the lack 
of correlations among the lengths lij brings an analytical solution within reach. 

It was noted by Vannimenus & Mezard (1984) that in the limit of large N, one may 
perfectly well choose the distribution of the individual Z^j's in the random link tsp so as to 
match that of the cZ-dimensional Euclidean TSP. The probability of finding a point at Euclidean 
distance I is simply the surface area of a d-dimensional sphere of radius I: 

In this way, the random link and Euclidean models can be made to match one another, up to 
correlations between lengths. We may therefore consider the random link TSP as a random link 



Alternatively, the final constraint can be expressed by considering a partition of the N points into two sets 
A and B, and requiring that V^, B, J^ieA jeB — ^■ 

■'Asyiiiiiiotric versions of the probloiri do exist as well, though we do not eoiisider tli(>m here. 
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approximation to the Euclidean tsp. This approximation is far from perfect, since in reahty 
correlations between Euclidean distances, such as the triangle inequality, are important. In terms 
of the results obtained, however, the random link approximation turns out to be a surprisingly 
good one. 

The approach in this thesis is, first of all, to examine the d-dimensional Euclidean TSP, 
focusing on the large N behavior of the optimal tour length Le ("e" for Euclidean). Our study 
is numerical, making use of good approximate methods for estimating Le to within small, well- 
controlled errors. This, along with some theoretical insight into the scaling behavior to expect 
in arbitrary dimension d, enables us to obtain (in units where volume equals one) the large N 
finite size scaling law for the ensemble average oi Le- 



(LE[N,i))=DE[d) N^-"" 



(1.4) 



At d = 2 and d = 3 we find extremely precise numerical values for /^^(d): our results, /3e(2) = 
0.7120 ± 0.0002 and /3b(3) = 0.6979 ± 0.0002, represent a precision 20 times greater than was 
possible using previously available methods.^ It would of course be nice to have an analytical 
estimate of /?£;(d), and this is where the random link approximation comes in. As wc have 
mentioned, this approximation is what makes the TSP analytically tractable. Using the "cavity" 
approach of Mezard, Parisi and Krauth (Mezard & Parisi, 19866; Krauth & Mezard, 1989; 
Krauth, 1989), we obtain an estimate for /3E{d). For d = 2 and d = 3, the random link 
approximation turns out to be a remarkably good one, with an error of less than 2%. In the 
d ^ 00 limit, we argue that in fact the approximation is exact up to 0(l/(i^). This leads to a 
conjecture on the behavior of /3b (d) at large d, in terms of a series in 1/cZ where leading and 
subleading coefficients are given. 

It is worth explaining why these Euclidean results are of broader interest than just to the- 
oretical physicists. One of the earliest mathematical results for the Euclidean stochastic tsp 
(Beardwood, Halton & Hammersley, 1959) states that the random variable Le obeys a 
self-averaging property: 

^ PE{d) as iV ^ 00 (1.5) 

with probability 1. /3Eid) is a well-defined quantity and not a random variable. This means 
that if we take a sequence of instances drawn randomly at large N, the relative fluctuations in 
the optimal tour length will be small from one instance to the next. (See Appendix B for an 
outline of the self-averaging proof.) Using this property, we may consider some large real- life 
Euclidean situation, such as an airplane serving numerous cities or a delivery truck serving a 
large number of retailers, and make the assumption that the instance "resembles" one from the 
random ensemble, i.e., sites positioned independently. Knowing PE{d) allows us at least a rough 
estimate of how far any prospective tour in this instance is from "theoretical optimality" . (A 
concrete example of this is given in Appendix C.) 



^These values have subsequently been confirmed by Johnson, McGeocii & Rothberg (1996), using a mod- 
ified version of our iiu^tliods and more powerful computational rosour(x;s. Sec Appendix E for further discussion. 
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From a theoretical point of view, the random hnk tsp is very much of interest as well, as 
a model of its own. The cavity method gives an analytical prediction for f^Riid), the random 
link counterpart of (3E{.d). However, this prediction is based on some rather subtle assumptions, 
reflecting work done in spin glasses. The central such assumption is that of replica symmetry, 
an ergodicity hypothesis stating that there is a unique thermodynamic limit, in a certain spin 
model onto which one maps the random link TSP. This assumption is known to be false in the 
case of spin glasses. For the tsp, however, numerical analyses performed in the d = \ case 
(SOURLAS, 1986; Krauth & Mezard, 1989) Suggest that it is in fact correct. We extend this 
numerical work to higher dimensions; our results indicate strongly that the cavity method gives 
correct predictions for all d. 

Finally, there is a related stochastic problem, discussed in this thesis, that does not directly 
involve the tsp though it is motivated by the finite size scaling analysis Tiscd for the ETiclidcan 
TSP. This is a geometric issue concerning properties of distances between randomly distributed 
sites, notably distances from a point to its A;th-nearest site. If the sites are distributed uniformly 
in Euclidean space (as cities are in our formulation of the Euclidean TSP), these fcth-neighbor 
distances display a remarkable universality: they all have the same large N scaling law. Letting 
{Dk{N)) be the ensemble average of distances between A;th-nearest neighbors, the A'^-dependence 
and the fc-dependence of {Dk{N)) separate, and 

[r(d/2 + l)]V'^ Tik + l/d) T{N) 

^'''^''^^ = — m~ m+i/d) ^'-'^ 

in d-dimensional Euclidean space. If we work not in Euclidean space but on an arbitrary 
geometric manifold, a different universality emerges: when {Di^{N)) is written as a power of N 
times a correction series in the 0{l/N) term of the series is a topological invariant and 

does not depend on the detailed properties of the manifold. We discuss these universalities and 
surrounding issues, noting that some open questions include the application of this work to more 
complicated geometric problems such as triangulations. 

• • • 

In both its stochastic and non-stochastic forms, the TSP has attracted a large amount of 
research. In order to provide the reader with a better perspective of the context in which our 
own work is situated, let us survey some of the most relevant results in the literature. 

The Utopian ideal for a TSP algorithmic procedure would be one which finds the exact 
optimum, for an arbitrary instance, in a short amount of time. The fact that the problem is 
NP-hard appears to make this goal unattainable (see Appendix A for a discussion of the TSP's 
complexity). Nevertheless, computer scientist have for years tried to develop algorithms that 
reduce the amount of time necessary to find a guaranteed optimal tour. The most successful 
such algorithms have been based on the method of "branch-and-bound" . The idea is as follows: 
a recursive procedure is used to break down an AT-city problem into smaller subproblems. Each 
subproblem is then "relaxed" so as to find a bound. These bounds are subsequently used 
to eliminate further subproblems from consideration. In this way, the algorithm successively 
"branches" divides problems into subproblems and ■"boimds" primes branches that 
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violate bound restrictions (Balas & Toth, 1984). Finding efficient branching and relaxation 
procedures are, of course, highly sophisticated problems in themselves. A variant on this method 
known as "branch-and-cut" represents today's state of the art in exact algorithms: problems of 
up to 100 cities can typically be solved in minutes on a workstation (Padberg & RiNALDi, 1987). 

However, when instances of larger sizes must be optimized, or a large sequence of instances 
must be optimized, exact algorithms are no longer practical. The next best choice, then, is to 
use approximate methods that find near-optimal tours quickly (typically in 0{N'^) operations). 
These are known as heuristic algorithms. For the TSP, the most powerful heuristics have been 
those that start with a non-optimal tour and successively attempt to improve upon that tour. 
The procedure generally used is known as local search: a certain number of links are removed 
from a tour and arc reconnected so as to form a new and better tour; when all possible local 
changes have been tried without success, the resulting tour is then a local optimum. The sim- 
plest algorithm of this sort is 2-opt, where two links are deleted and the paths are reconnected 
differently. 3-opt does the same but using 3 links in each elementary move (LiN, 1965). Both 
heuristics tend to give tour lengths that are a few percent higher than the optimum for Euclidean 
TSP instances (JOHNSON &; McGeoch, 1997). Various methods have been popular for improv- 
ing on this. The simplest is to repeat 2-opt or 3-opt from diff'crcnt random starting tours, taking 
the best local optimum found over many trials, so as to cover more of the energy landscape (see 
Figure 1-2). Another method is simulated annealing, which can outperform repeated 3-opt (in 
terms of both running time and percent excess of optimum) only for sufficiently large instance 
sizes {N > 1000). The most successful of the local search heuristics, however, is that of LiN & 
Kernighan (1973). This algorithm uses a "variable-depth neighborhood search" method where 
a large number of links can be changed at once, but where redundancy in searching is reduced 
to the point where a million-city Euclidean instance can be solved in less than an hour, with 
less than 2% excess over the optimum (see Appendix D for a description of the heuristic; see 
also Johnson & McGeoch, 1997). 

Improving upon Lin-Kernighan has been a goal in algorithmic research for over 20 years, with 
relatively limited success. One particularly creative attempt has been through the use of genetic 
algorithms where, inspired by the process of evolution, one attcmps to establish new starting 
tours by "mating" (combining aspects of two "parent" locally optimal tours) and by "mutation" 
(randomly altering a locally optimal tour). Although the mating principle has not yet proven 
to produce a competitive heuristic, the mutation principle has been used with success in the 
"metaheuristic" known as Chained Local Optimization (cLo) by Martin, Otto & Felten 
(1991). CLO involves successively applying a random "kick" to a locally optimal tour, and then 
using Lin-Kcrnighan to locally optimize from there. If the new tour is better than the old one, 
the procedure iterates from that point; if not, the old tour is again used and a new random 
"kick" is attempted. (A qualitative description may be found in Appendix D.) CLO outperforms 
Lin-Kernighan for instances of iV > 100, taking seconds to get within 1% of the optimum for 
a 1000-city instance as opposed to minutes for the equivalent performance simply repeating 
Lin-Kernighan from random starts (Martin & Otto, 1996). 

All of what has been mentioned so far applies to optimizing individual TSP instances. When 
we speak of the stochastic TSP. wo must access an entire ensemble of instances. One approach 
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is simply to run optimization heuristics over a large sample from the ensemble, thus obtaining 
results to within a quantifiable statistical error. This approach is indeed the basis of most numer- 
ical studies, including our own. These numerical studies, however, have also been accompanied 
by a certain amount of theoretical work. 

For the Euclidean stochastic TSP, we have already seen the self-averaging property shown 
by Beardwood, et al. These authors, in the same article, gave lower and upper bounds on the 
asymptotic limit PE{d)- The lower bound is simple: it uses the fact that, at best, any city on 
the tour has links to its two nearest neighbors. Nearest neighbor distances may be calculated 
over the ensemble (as we discuss in depth in Chapter IV), and one finds that /3e(2) > 5/8. It 
is straightforward to generalize this bound to any dimension d. For the upper bound, a method 
known as "strip" was used, whereby the space is divided into columns (strips) and a suboptimal 
tour is constructed by visiting cities alternately up and down the columns. In 2 dimensions, this 
suboptimal tour length may be calculated exactly over the ensemble; the result is a complicated 
expression involving special functions, but may be numerically computed as /3£;(2) < 0.9204. 

No better exact bounds are known today; the best that has been done is to use heuristic 
methods on very large instances (or many very large instances) to establish estimated bounds. 
In this way, an approximation to the lower bound of Held & Karp (1970), involving a linear 
programming relaxation of the TSP, has very recently been used by JOHNSON, McGeoch & 
ROTHBERG (1996) to estimate PeC^) > 0.708. Lee & Choi (1994) have used a variant on simu- 
lated annealing to obtain the upper bound /3£:(2) < 0.721. At c? > 2, virtually no improvements 
over the exact bounds have existed prior to our own work. (Johnson et al. have, however, in 
their work on the Held-Karp bound, successfully applied a modified version of our numerical 
methods to higher precision, confirming our own d = 2 and d = 3 results and extending to the 
d = 4 case.) The only other related large d result is an unproven conjecture by Bertsimas & 
VAN Ryzin (1990), claiming that when d — > 00, PE{d) ~ y^d/2TTe. 

For the random link stochastic TSP, results have been largely due to the work of theoretical 
physicists. Vannimenus & Mezard (1984) and KiRKPATRicK & Toulouse (1985) first con- 
sidered the statistical mechanics of the TSP, using analytical tools stemming from the study of 
spin glasses. Vannimenus and Mezard showed that the equivalent of the Bertsimas-van Ryzin 
conjecture is correct for the random link case, i.e., PRL{d) ~ \/ d/2TTe when d 00. Kirkpatrick 
and Toulouse studied several low-temperature macroscopic quantities of the random link model, 
analyzing the number of links in common between local optima at finite temperature. This 
"overlap" leads to the notion of a state-space bond distance, measuring by how much two tours 
differ from each other. In analogous spin glass problems, such distances between local minima 
display a property known as ultrametricity , a stronger form of the triangle inequality, requiring 
that for any three distances at least two must be equal (they form an isocelcs triangle) and the 
third must be less than or equal to the other two. Non-trivial ultrametricity is considered a 
typical sign of a lack of ergodicity in spin glasses, called replica symmetry breaking. SOURLAS 
(1986), numerically extending Kirkpatrick and Toulouse's low-temperature analysis, found on 
the basis of the tour overlap that ultrametricity does not play a role in the d = 1 random link 
TSP, and that replica symmetry most likely holds. In the process, he observed that solutions 
only contain links between very near neighbors, and used this property to motivate an improved 
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optimization heuristic that works by disaUowing certain hnks. 

Under the assumption that the rephca symmetric solution is the correct one, Mezard &; 
Parisi (19866) and Krauth & Mezard (1989) used the cavity method to find an analytical 
expression for f3RL{d), in the form of an integral equation. They performed further numerical 
checks at d = 1 (through direct simulations of the random link TSP) that confirmed their 
analytical cavity results, and briefly discussed the d = 2 result. We are aware of no further 
attempts, prior to our own, to extend the analysis of the cavity solutions to higher dimensions. 



Let us conclude this chapter with a brief guide to the layout of what follows in the thesis. 
Chapter II deals with scaling laws in the stochastic TSP. In the articles presented in the chapter 
we discuss our most important results, notably a numerical study of the finite size scaling of 
Le^ and an analytical study of the dimensional scaling of Psid) by way of the random link 
approximation. Our analysis of LE{N,d) at large is, to our knowledge, the first to use 
a systematic method for extrapolating to the large N limit. We are thus able to improve 
significantly on previous estimates of fiE{d). Using the random link approximation, we also 
present evidence supporting the Bertsimas-van Ryzin conjecture that PE{d) ~ ^/d/2ire at large 
d, and propose a much stronger conjecture: that 



where 7 is Euler's constant (7 0.57722). 

Chapter III covers the random link TSP. We discuss the background of the analytical tools 
used, and the physical model serving as the basis of the cavity method. Through simulations, 
we then test the cavity results for d > 1, and find good numerical justification for the hypothesis 
of replica symmetry. 

Chapter IV provides a discussion of the universality in fcth-nearest neighbor distances for 
randomly placed sites in Euclidean space, namely that the same scaling in N applies to the mean 
distance {Dk(N)) regardless of k. We generalize this study to arbitrary geometric manifolds, 
and find that, although the scaling of {D^lN)) does depend on k in general, its 0{1/N) term 
is a topological invariant and does not depend on the precise shape of the manifold. We note 
that these properties, though arising from relatively simple physical arguments, can be applied 
to far more complex geometric problems. 

In each of these subsequent three chapters, an overview of the subject is given, following 
which one or more articles make up the body of the chapter. In Chapter II, these are articles 
that have already been published, and are presented in reprint form. In Chapters III and IV, 
the articles have not yet been submitted, and so it has been possible to take some liberties with 
their presentation, following wherever possible the style adopted in the rest of the text. Not 
being constrained by editorial restrictions concerning brevity, we have also attempted to adopt 
a more pedagogical tone than will ultimately be permitted in published form. 

Several appendices follow the main body of the thesis. In Appendix A, we define the notions 
P and NP as used in computational complexity theory. Wo show how to phrase the TSP in such 
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a way that it is an NP-hard problem, and discuss polynomial time approximation algorithms. 
Appendix B covers self-averaging in LE{N,d); we outline a simple proof of this property, es- 
tablishing that LE{N,d)/N^'~^^'^ goes to a constant ^Eid) when N ^ oo. In Appendix C we 
explain how our Euclidean stochastic TSP results, although intended for a uniform distribution 
of cities, may easily be applied to non-uniform cases. We give the example of the AT&T-532 
instance (containing 532 AT&T telecommunications sites around the U.S.), and show that we 
come within 1% of the true optimum by regarding it as a member of the stochastic ensemble. 
In Appendix D, we discuss our numerical methodology, giving an explanation of the workings 
of the Lin-Kernighan and CLO algorithms for optimizing TSP instances. Finally, in Appendix E, 
we describe a more recent numerical study of (3Eid) and (3jiL{d) performed by JOHNSON, Mc- 
Geoch & ROTHBERG (igg6), using a variant on our methods. They obtain results in close 
agreement with ours, and by combining their data and ours we are in fact able to refine our 
PeC^) estimate slightly further. 
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— Chapter overview — 



Historically, the traveling salesman problem has been of greatest interest in its Euclidean 
formulation. In this form of the problem, each of the N cities has a position in d-dimensional 
space (in practice, most often 2-dimensional space), and the lengths separating pairs of cities 
are defined to be the associated Euclidean distances. The traveling salesman tour is thus a tour 
through a metric space. 

Much past work on this problem has concentrated on developing algorithms to solve certain 
large instances, or often certain classes of instances. One goal in operations research has been 
to find methods of reducing redundancy in exhaustive search — methods that, while still expo- 
nentially slow in N, at least have a small enough coefficient in the exponential that they find the 
exact optimum in a reasonable amount of time. Another goal has been to classify special cases 
of instances. The idea here is to identify properties of the layout of cities that could guarantee 
finding an optimal tour within, say 0(N"') steps for some fixed a. Finally, considerable efforts 
have been devoted to finding mathematically rigorous bounds on tour lengths. 

In the case of the stochastic Euclidean TSP, however, the emphasis is somewhat different. 
The positions of the cities in the tour are now independent random variables; in the model that 
we consider here, the distribution is uniform over space. ^ It is no longer a matter of optimizing 
a single instance. We are now working with an entire ensemble of instances, and the problem 
is to find properties of the optimal tour length Le (the actual path that tours follow is of less 
interest when the positions of the cities are not fixed quantities) considered in the ensemble of all 
possible instances. The challenge consists of understanding what happens at large instance sizes. 
The fundamental property here is that of self-averaging , proven by Beardwood, Halton &; 
Hammersley (1959): as N gets large the relative fluctuations in the optimal tour length go to 
zero, and the distribution of Lg, up to a scaling factor, becomes more and more sharply peaked. 

^The non- uniform case is discussed in Appendix C. 



Chap. II: Scaling laws in the Euclidean TSP 



The real quantity of interest is then the mean optimal tour length over the ensemble — at a 
given instance size N and for a given dimension d — which we denote {LE{N,d)), using units 
where the volume of our space is 1. 

An outline of a more recent and simpler proof of self-averaging, due to Karp & Steele 
(1985), is given in Appendix B. The precise statement of the property is that with probability 1, 
for any sequence of random instances, Le{N, d)/N^~^^'^ will approach the instance-independent 
quantity Psid) at large A''. Note that the scaling factor N^~^/'^ may also be seen by a simple 
physical argument: we are working in a fixed volume, so the mean volume per city scales as 
N"^, and so the mean distance between cities scales as N~^/'^; the tour contains N such links, 
hence its length scales as N^~^/'^. 

The question then becomes how to determine the limiting 

For small values of d, this may be answered by combining numerical methods with some theoret- 
ical insights — namely, how LE{N,d)/N^~^/'^ converges to its large N limit. We thus consider 
the finite size scaling of LE{N,d). At large d, this no longer feasible numerically. Fortunately, 
in the limit d 00, an analytical approach may be developed for obtaining f3E{d) to within an 
excellent approximation. 

What sort of numerical methods are needed for a finite size scaling analysis? At a given 
d, we must find {LE{N,d)) over many values of N. At any particular value of N, in order to 
obtain this quantity we must, in turn, average Le{N, d) over a large number of instances. (The 
smaller N is, the more instances we will have to consider in order to minimize the effects of 
fluctuations.) At any given instance, finally, in order to find the optimum we must have efficient 
algorithms. Unfortunately, exact algorithms are too slow to be of use here, requiring minutes of 
CPU time to solve typical instances of N = 100; optimizing a statistically significant sample of 
random instances in this way could take weeks of computing time. 

There is, fortunately, a better approach. Powerful heuristic algorithms exist — inexact 
algorithms which find "good" tours but not always the optimal one. The most suitable method 
is that of local search, where the algorithm takes a non-optimal tour and iteratively improves 
upon it so as to bring it closer to optimality. We may think of this in terms of the landscape 
in Figure 1-2; the algorithm finds the local minimum of the valley in which we start, and if it 
is a good algorithm, the landscape will be such that there are a small number of broad valleys. 
As we execute a local search heuristic increasingly many times on a given instance, each time 
starting from a different random initial tour, we then explore increasing parts of the landscape 
and find the true optimum with increasing probability. Moreover, even with multiple runs, 
these heuristics work considerably faster than exact algorithms — especially at large A'^ where 
running time grows only polynomially and not exponentially in N. This means that in a given 
amount of time, we are able to optimize a much larger sample of instances (and thus have 
a far smaller statistical error) using heuristics than using exact algorithms. The net effect is 
counterintuitive: because of computing time considerations, we actually get more precise results 
using these inexact algorithms than Tising exact algorithms. There is. admittedly, a systematic 
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bias in the results arising from the use of quasi-optimal tour lengths: no matter how many times 
we run a heuristic, we are never sure to have found the global optimum in the end. However, 
we choose the number of runs per instance sufficiently conservatively that this systematic bias 
— estimated for each value of A'^ on the basis of a number of test instances — is so small as to 
be negligible compared with the statistical error. (Our numerical methodology is discussed in 
Appendix D.) 

Equipped with this numerical approach and certain theoretical notions concerning the scaling 
to expect in N — described in detail in the articles that follow — we find the large N finite size 
scaling law: 



using a d-dimensional unit hypercube with periodic boundary conditions.^ Understanding the 
scaling in this manner enables us to obtain Pe (d) to high precision at d = 2 and d = 3 (the cases 
that we consider numerically): /3e(2) = 0.7120 ± 0.0002 and /3e(3) = 0.6979 ± 0.0002. When d 
gets larger, however, the feasibility of numerical solutions rapidly decreases. In order to get an 
idea of the dimensional scaling of /3b (d), we therefore turn to the analytical approach used in 
the random link approximation. 

The random link approximation consists of assuming that lengths lij between city i and city 
j {i < j) are completely independent of one another. Thus under the random link approxi- 
mation (which could perhaps more properly be called the "independent link approximation") 
correlations between distances, such as the triangle inequality, are neglected. Making this ap- 
proximation allows us to take advantage of an analytical approach developed by Mezard & 
Parisi (19866), Krauth & Mezard (1989) and Krauth (1989). We use this approach to 
obtain fim {d) , the random link value approximating Pe (d) ■ Prl (d) arises from a system of in- 
tegral equations and does not appear to have an exact analytical solution, but can on the other 
hand be solved numerically to arbitrary precision. For small values of d, we find that Prl (d) is 
a good approximation to f5E{d): for d < S, the relative error is within about 2%. 

For large d. we are in fact able to express (3rl (d) analytically, via a power series in 1/d. To 
leading and subleading order, this gives 



where 7 represents Euler's constant (7 0.57722). We then present an argument, based on 
a theoretical analysis of how certain heuristic algorithms operate, suggesting that the relative 
error in the random link approximation itself decreases at least as fast as 0(l/d^) in the limit 
d — > 00. This enables us to conjecture, finally, that the large d expression (II. 3) applies to PE{d) 
as well, and thus gives the correct dimensional scaling for the Euclidean optimal tour length. 



^Jaillet (1993) has proven that (3E{d) is in fact the same regardless of whether periodic or open boundary 

conditions arc used. 




(11.2) 
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Finite Size and Dimensional Etependencein the Euclidean Traveling Salesman Problem 

A lion G. Be Lieu s* and Oliviei' C. Mail in* 

Drviiioii de fhfi\qss Theouqae, hiimse ds fkfikqae Nucisairs, U)irveriils Parii-Sud. F-914G6 Ormy Cedsx. Fra)ice 

tRccciird^O Jijl> 1995) 

Wc oonEidjCL the Euclidean tfa deling EalcEman pLToblcm foi N citicE landjomlj distnbiited in the 
unit ^-dimcnEio nal hjpcL"cube, ind in^xEtigite the finite size Bcahng of the mean optimil tout" 
length Lf. With toroidal boundaij conditior^ we find, mofivated by a lejn aftab le univeiEalitjf in 
the Ath neaiest neighboi" diriii bution, that Lcld = 7) = (0.7 IK) ± O.OOffJ) W '"^ [ I + Oll/N)] and 
LfW = 3) = (0.6979 * 0.00O2)W^^I + 0[l/N)]. We then coraidei' a mean-field appi-oach in the 
limit N ^ ^ which we find to be a ^od appicximafion [the eiroi" being Icee than 2.1% at ^ = 1,2, 
and 3), and which Eug^EtE that Lf(<f) = N' ^^''.Jd/lur (-jTifl'-^Cl + Oll/d}] at laige d. 



FACS numteis: 02.60.Pn, 02.710.Ljq, MjMjCn 

The tiareling saljesLnan piDbienn (TSP) is one of the 
best known corabinalDLial optimisation problems. It is 
NP complete (suggesting that no algoiithm eiists fof 
solving the pioblem in polynomial time), and it sen/es as 
a feitile ground fof analytical and numeiical appiioaches to 
optimization problems in geneiai. It is also one of the few 
optimization pioblEins tha hai/e been studied eitensively 
in the oonteit of statistical mechanics. 

The TSP, as we con side f it, is as follows: Given A/ 
points ("cities") in a space, the problem is to find the 
length of the shortest closed path ("tout") going thtough 
each city exactly once. Two particulai' fbi ms of the prob- 
lem have been investigated in depth. The fiist, which has 
attiacted the most attention among computet' scientists and 
nnathennaticians, is the Euclidean TSP: The N cities ate 
tandomly disttibuted in a i^-dinnensional hypeicube and 
the distances between cities aiie given by the Euclidean 
met lie. The second, which has been of particulaf inteiiest 
within the statistical physics connmunily, is the landonn 
link TSP: The lengths fjj sepaiating cities / and J ate 
taken as independent landonn vatiables with a given dis- 
ttibution pit). 

Tt has been noted by M^zaid and Paiisi [1] that the 
tandom link model, with pi!) approptialely chosen, maps 
onto the Euclidean model if cotiielalions between thiiee ot 
nnoLie distances aiie neglected (no ttiangle inequality, fot 
instance). This suggests that the landonn link TSP can be 
con side Lied as a mean-fieldappiDiimationto the Euclidean 
case, and peihaps that this sppioiimation becomes eiact 
in the liiTiit d 

Out intention in this Lfittet is twofold. Fitst, fot 
the Euclidean TSP we investigate finite size cotiections 
to the mean optimal tout length , in the laige N 
("theimodynamic") limit. To out knowledge the lie has 
been no ptiot woik on this subject, in spite of a giieat 
deal of inteiiest in in the thetnnodynannic linn it itself 
Second, we eiploiie the dinnensional dependence of 
using a mean -fie Id appiioach (the tandoiin link TSP in 
conjunction with the "cavity method" [U'J]). We extend 
the woik of Ktauth and M^zaid [3] to find the nnean- 
field optimunn L^jp in the thetmodynannic I inn it, as a 

1188 003 1 -90af7/96/76(8)/ 1 1 88(4)i06.00 



function of dimension. Compaiing mean -fie Id tesu Its with 
Euclidean N lie suits at low d shows that mean field 
does consideiably bettet than piieviously expected, and 
suggests that in quite natuial units, can be w title n as a 
powet seties in I /d. 

Euclidean mockl: Finite si^ scaling Id = 2]. — We 
stait with the case of A/ cities disttibuted tandomly 
and unifotmly in a unit squaiie. Nunnetous heutistic 
appiioaches have been developed to find neat-optimal 
TSP touts given a particulat configutalion ("instance") of 
cities. Fot out puiposes, the most convenient methods ate 
local-optimization heuiistics such as the Lin-Keinighan 
(LK) [4] and the chained local optimization (CLO) [5] 
algotithms. With these algptithms, lepeated mns on a 
given instance using diffetent tandom starts produce the 
optimal tout with inc teasing ptobability. 

Tt has been shown [6] that in the laige N limit the 
optimal touf length jbr a gitvn institnce is self- 
aveiaging up to a scaling factot 

wheiie conveigence to the instance -independent is with 
ptobability 1 (in the ensemble of instances with tandomly 
disttibuted cities). Much past woik has concenttaled on 
optimizing single instances at laigp N [see [5,7,8]). He lie, 
ho we vet, out concern is to calculate along with an 
estimaiE of statistical etrot, and so Instead we aveiage 
ovef a laige numbet of instance. The lie is necessaiily 
a tiadeoff in the choice of A/: At small A/ alone we 
cannot confidently ptedict the finite size scaling behaviot, 
wheiieas at laige N the laige annount of computing time 
necessaiy fot each optimization shaiply limits the numbet 
of instances we can optimize teliably, and incieases the 
statistical ettot. We thetefote choose seveial small values 
of A/ (A/ = 1 5 through A/ = 17) w he iie we optimize using 
LK, and two laige t values (N = 30 and N = 100) wheie 
we optimize using CLO. 

Given (A/) at diffetent values of A/, then, we wish 
to extiapolale and extiact the limit , as well as finite 
size collections. In oidet to eliminate the effects of 

© 1996 The Ameiican Physical Society 
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suiface teiiTis, we use pei iodic boundaiy oonditions in 
the Euclidean distance metiic. An indication of the size 
dependence to be eipected in L^IN) wa^ be found 
by looking at the distance Dt between klh neaiiest 
neighboLS, ai/Ei'^ed ovei' the enseiuble of instances. A 
diiect calculations shows that, given W citijes distiibuted 
landoiTily and unifbiiiily oi/ei' the ^^-dilTensional unit 
hypeicube (with pei' iodic boundaiy conditions). 



IN - k)d 



(:::) 



+ I) 

IT-'/ 



] 



wheiie eiponentially small collections in N hare been 
neglected. 

Recognizing this integral (up to a change of vaiiable 
and funhei' eiponentially small collections in A/) as a beta 
fiinction, we find that 



DilN,d] 



TIN + \/d] Vif 



(1) 

Notice that theie is a complete sepaialion heie of the 
A/ dependence and the k dependence. This is indeed 
a suipiising unii/ieisality: Tt means that up Id eiponen- 
tially small collections, all kth nearest neighbor mean 
distances hare exictly the same scaling law in N, namely, 
riN)/riN + \/d). It might be eipecled, then, that the 
length of a TSP toui' consisting of N links would have 
laige A/ scaling behavioi' 



N 



T{N + \/d] 



I - yd 



I + 



i/d - i/d- 



wheiie the light- hand side follows fiiom Stiiling'sfoimula 
Tn fact, due to coiiielalions between k and hi in the 
optimal toui', this is not quite the case. Figuiie 1 shows 
oui' iiesultsfoi' divided by the scaling quantity above, 
sld = 1: We find that this is, to a good fit, itself a powei' 
seiies in \/N, albeit one with a small fiist-oidei' teim. 
The asymptotic N ->■ value is = 0.71-20 ± 0.0002, 
wheiie the eiioi' is obtained on the basis of x' analysis. 
This lesult is, to oui' knowledge, the most piecise to dale 
fbi' the Euclidean TSP in the theimodynamic limit. 

The methods by which we obtained the lesults in Fig. 1 
aiB themselves of some importance. Foi' luns optimized 
by LK {hi = \1 thiough N = \7), we aveiaged ovei" 
the iiesults of ^50 000 instances, wheie foi' each instance 
we took the best (lowest) optimum found in ten landom 
starts (ten different nins). Foi' A/ = 30 we aveiaged ovei' 
10 000 instances, taking fbi' each one the best optimum 
found by CLO (ten Monte Cailo iteiations pei' mn) in 
five random starts. Foi' N = 100 we aveiaged ovei' 
6000 instances, taking fbi' each one the best optimum 
found by CLO (ten Monte Caiio iteiations pei' ran) 



0,7126 



0.7105 



^ 0.7085 



0.7065 - 



0-70+5 



0,7025 




1/W 

FIG. 1. Finite size dependence of iieECiied Euclidean 2D 
TSP optimuiri. Best fit = SAS) is given \sy Lc/N'^-ll + 
1/[SN) + ...] =0.71-50(1 -0.0I7I/JV - imS/N^. Eiroi- 



in twenty landora starts. These methods intioduce a 
systematic ernoi', because they do not always find the 
tiue optimum; we estimated this eiioi' by peifbiming a 
laige numbei' of nins on a few instances and raeasuiing 
the avei'^e eipected eiioi' [weighted by the probability 
of making that eiioi' when choosing the best out of 
ten landom staits). In all cases, we veiified that the 
systematic eiioi' stayed undei' 10% of the statistical eiioi' 
shown in the eiioi" bais. 

la oidei' to iieduce the statistical noise furthei', we used 
the following vaiiance iieduction method: Recognizing 
that Lglhl] = NID, + Di)!! is a bwei- bound on the 
toui' length (each city is at best connected to its fiist- 
and second -neaiiest neighbois), wiite the estimatoi' fbi' 
as {Le - ^Lg) + XLg. Le and Lg denotes values /or 
a particular instance, the angulai' biackets iiepiesent the 
aveiage ovei' instances saiuple, and the ensemble average 
Lg can be calculated analytically [see Eq. (1)]. A is a 
paiametei' which we adjust to minimize the vaiiance of 
oui" new estimatoi". Tn piactice, optiitial values of A (A 
0.75) enabled us to iieduce the ei ioi' by ovei' 60%. Othei' 
vaiiance reduction methods can also be used [9], but ouis 
has the advantage of intioducing no new systematic ernoL 

yfean-Jield ttiethad. — We now tuin oui' attention to the 
mean-field appioiimation, based on the landom link TSP. 
Rathei' than having N cities distiibuted landomly in a 
hypeiicube, we now have lengths fij between cities / and 
j (I ^ i < j hi] distiibuted as independent landom 
vaiiable s accoiding to a certain distiibution fill). We 
take fill) to be the probability distiibution of lengths 
between cities in the -dimensional Euclidean pioblem, 
in the absence of finite size effects: 



fill] = dir''f-l''-'/rid/2 + I) 
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This establishes a mapping in the theiraodynaraic limit 
between the Landora link TSP and the Euclidean TSP, 
neglecting all coiielations among (Euclidean) distances. 

The LTiean-field "model" is the landom link TSP, 
desciibed fot oui' puiposes by the "carity equations" 
wi itten down by Kiauth and M&and [3]. In oui" langu^e 
this leads to 



rjj + I) 

wheiie /^uF Lh^f/^'" In the Euclidean case, and 



It has been aigued peisuasively, notably on the basis of 
excellent agreement in the = I case [3], that the cavity 
method is enact fof the N — * lanidom link TSP. In 
the following discussion we shall also piiesent fuithef 
justification foi' this assuiTiption. 

Theie is no known analytical solution of the integial 
equation foi" GAx] given in Eq. (^I Howevtei", it can be 
solved numeiically; this was done by Kiauth and M^zaid 
at = I and d =2, giving ^n^ld = I) = I.020S and 
^nr^ld = 1] = 0.725 i [3]. These values may be cora- 
paied with field): Undef peiiodic boundaiy conditions 
fisld = \ ] = i [tiivially) and fi^id = 2) = Q.li2f> (see 
piievious section). Theiefoie, ax d = i mean field has a 
2.1% excess with respect to the Euclidean value, and at 
d = 2 a 1.8% excess (see also Table I). Aliieady at low 
dimension, then, mean field gives quite a good sppioxi- 
matbn to the Euclidean case. It is amusing to note that 
Kiaut handM^aidthemselvesassumedaiatheii naccuiale 
Euclidean value field = 2] = 0.749, and so theii' mean- 
field results seemed poorei lotheiT than they actually weie. 

We now extend the numeiical solution of Eq. (2) to 
highei dimensions. As in the piioblem of Euclidean 
finite size scaling, we can get an indication of what 
dimensional dependence to expect in Lnf:ld] by looking 
at the mean kth nearest- neighbot distance Di multiplied 
by the numbet of links A/. In the theiraodynamic limit, 
Eq. [1 ) gives 



NDild) 
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i/d) 



m) 

d 

f 2lTf 



at Laige d . 



Dividing by this suggests that 



fild) =J—i-!Td) 



I + 



TABLE I. CojnpiLiEon of Euclidean and mtan-field TSP 
optinn-ii fL"CECiilcd) -at diincnEion up to -tf = 3. 



d 








1 


1 


1.020S 


+ 2 l?S. 


2 


0.7 110 ± 0.00012 


0.725 1 


+ I.B?5i 


3 


0.6979 ± 0.0002 


0.7 100 


+ 1.7% 



Figure 2 shows that this is indeed so fbf the mean- 
field results obtain ed by n umei ical lesolution of Eq. (2). 
Looking at fifjif/-Jd/2iTe ( ird)'^-'', we find an excellent 
fit by a \/d powet seiies with a leading oidet teim 
which, to the piecision of ouf law numeiical data, is 
indistinguishable ficim 1. 

The faa that finr'/-.Jd/2ire at d ^- is anoihei- 
confiiraation of the validity of the cavity method, as this 
piopeity is known to be tme foi the puie landom link 
TSP [10]. We have thus added to Kiauth and M^zaid's 
investigation (at = I ) furthei evidence (at li — » a^^) that 
the cavity method is exact. 

Finally, let us new lite the left-hand side of the best-fit 
equation in Fig. 2 with an additional l\/2]'^^ facloi in 
the de nominator 



fin^ 



-Jd/2ireliTd/2] ^i-" 



= 0.999997 + 



0499395 



+ 



<7) 



Notice that the I /d coefficient is piactically indistinguish- 
able ftom \/2. An inteip fetation of this re maikable result 
is given in [1 1]. 
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2.4 



2.2 




FIG. 3. RcEcaled Euclidean TSP optimum fpoints) as a fij Fic- 
tion erf dimension, Eindwiched between mean-field optimum 
(solid line) and exact lowei' bound [dashed line). 



Eiiciiikan tttoslei: Ditnensiomi depenslen/x. — Gii/eri 
the LTieari-field lie suits, we now letuin to the Euclidean 
model. Table 1 shows the nu me ileal lesult at d = 3 
(obtained by the saiTie beuiisilc methods as in the d = 1 
case) tpgethei' with the mean -fie Id value, and the (J = I 
and = 2 IB suits piesented eailiei: 

These lesults suggest that §f4f Is an uppei bound foi' 
[and heuiistlc aiguments [11] also piciTlde support foi' 
this). At the same time, thei"e Is the stiict lowei" bound 
Lg = N{D\ + £>:)/"2, mentioned eaillei' In the discussion 
on vaiiance induction. Flguie 3 shows the Euclidean le- 
sults "sandwiched" between the coiiiespondlng mean field 
and lowEi bound quantities, both o f which may be wiit- 
ten In xVifid-^^ limit as fiid) = -Jd/lire liid] '■'^[1 t 
0(1 /d]'\. We conjectuiie that mean field does Indeed le- 
main an uppei' bound at all values c£d, and consequently 
that behaves asymptotically at laige d as 



^Eid] 



I + 



This would also support awealiEi' conjectuiie by Bert si mas 
and van Ryzin [12], staling that foi' the Euclidean TSP, 
— -.Jd/lire St d 
Jn conclusion, we fiave Investigated tfie finite size 
behavloi' of the Euclidean TSP optimum undiet peilodlc 
boundaiy conditions, and have seen that at = 2, 
conveiges as a I/A/ seiies: 

Le ^ ! . 0.0171 



Tn the piocess we have eiti' acted what we believe to 
be the best result to date fbt the theimodynamlc limit: 
^Eid =2] = 0.71201 0.0002. 

At the same time we have, by means of a mean-field 
method, eiamlned the dimensional dependence of the TSP. 
We have found that mean field Is a good appioxlmation (< 
1. 1 % eiiDi) to the Euclidean TSP m.d = 1 , 2, and 3. We 
have seen numeilcally thataliJ— '^i^the cavity equations 
aiie compatible with the eiact landom link TSP lesult, 
and thus have piDvlded fiirthet evidence that they aie 
eiiact at all dimensions. Additional woik Is In piogiess 
to undeistand the coefficient 1/2 in the subieadlng teim 
of the cavity equation solution. Finally, compaiing oui' 
mean-field and Euclidean results suggests not only that the 
Beitsl mas-van Ryzln conjecture fbt the laige d limit of 
^Eid) Is coirect, but also that the asymptotic behavioi' Is 
lnfact^£(ii) = -Jd/2se{'sd)^^-''[\ + OiXfd)']. 
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reseaich. We also acknowledge fmltfiii discussions with 
E. Bogomoiny, M. M&and, S. Otlo, and N. Souilas. 
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Note added. -S\.nct oui submission, D. S. Johnson nf. 
[13], using slightly dlffeient methods, have found values 
fbi' fisid) compatible with ouis at = 2 and 3. 
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The Random Link Approidniation for the Euclidean Itaveling 
Salesman Frohlem 
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(Rsasivsd 22 April 1&96, rarissd 2S Juna s.caaptad 2 A Saptamber 

PACIS.Q2.60.PiL Numerical aptimixatian 
PACf;.Cl2.7l].Lq Monte Carlo and atatistical nwtboda 

PACf?.64.6(].Cn Order-diaorder tranaldrmationa; atatdatical mechanica of model ayatema 

Abat lEnt . The traveling aaleamon problem (Tf!P) oonaista of finding the length of the short- 
est dosed tour visiting N "dtie^ . We consider the Eludidean TSP where the dtiea are dis- 
tributed randomly and independently in a<j-dimenaional unit hypercube. Working with periodic 
boundary oonditiona and inspired by a remarkable univeraslity in the feth nearest nei^bor dis- 
tribution, we find for the averaj^B optimum tour len^h (Lk) = .fls (d)JV'-'-^ [1+ 0(1/JV)] with 
,Sjc(a) = 0.7130 =t O.OOOa and ,3ic(,1) = 0.eS79 ±0.0003. We then derive analytical predictions for 
these quantities usin^the random link appraximatian, where the lengths between dties are taken 
as independent random 'wriables. From the "cavity" equations developed by Krauth, Mdicard 
and ParisL, we calculate the associated random link values ,S^{dy For d = 1,3,.1, numerical 
results show that the random link approximation is a ^od one, with a discrepancy of less than 
2.1% between ,Si£(d\ and ,Sai_{d). Fbr hrgo d, we ar^e that the approximsiion is exact up 
to Oll/<^) and ^ve a conjecture for jSatd), in terms of a power series in 1/d, specifying both 
leading and subleadin^ ooeffidents. 

nSauuiiS. Le probleme du vqya^eur de commerce (TSP) consisise a trouver le diemin fermi^ 
le plus court qui relie N "ville^ . Nous Studious le TSP eudidien ou les villes sont distributes 
au haaard de maniere dteorr^lte dana I'hypercube de cot^ 1, en dimenaion d. Fn imposant dea 
fxinditiona aux borda p4iriodiquea et ^id&a par uiieuiiiveraalit4 remarquable de la distribution dea 
feiemes voisins, nous trouvons la longueur moyenne du chemin optimal (L^) = JV^-^-'"' [1 + 

0(1/JV)], avoc ,9^(3) = 0,7120 ± 0j0002 et ,9ic(,1) = 0,6979 ± 0j0002. Nous fetablisaons ensuite 
des predictions analytiques sur ces quantitfas a I'aide de 1 'appraximation de liens al^atdres, ou 
laa lon^eura entre lea vdllaa aont dea mriablea altetoirea ind^pendantaa. Crace aux 4quationa 
"cavit^^ d4ivelopptea par Krauth, MAzard at Pariai, noua obtenona dana le caa de Liana alfaatoiraa 
laa valeura, analo^aa a ,8i£(d). Pour d = 1,2,.1, laa r^ultala num^riquaa oonSrment que 

Tapproximation de liens al^atcnres est bonne, conduisant a un 4cart inMrieur a 2,1% entre .Smld) 
et ,SiitL{d). Pour d ^and, nous donnons des ar,^uments montrant que cetlse approximation est 
exacte juaqu^a I'oidpB 1/^ et noua propoaona une conjecture pour ,8^(d), exprimte en Ibnction 
d'une a4rie an l/d, dont on donna laa deux premiera ordrea. 
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1 la JOTIHN AL DE PHYSIQUE I N°l 

1. Intraduction 

Giwn N "citira'" juid the diatJULcra between theiii^ the txavdin^ sAlefiiiiaii problem (TSP) 
eonadfitsi of finding the lengthof the sihiorteiit doaed "tonr" (path) visiiting every dty exactly on*x, 
where the tonr length k the sma. at the dty-tc-city difitaiu:ea pdang the tour. The TSP k NP- 
compktc, which anggeatsi that there la no general jilgorithiii capjible cf finding the cptimnni tonr 
in An Amonnt of time polynomial in N. The problem is thm aimple to abite^ bnt veiy difficnlt 
to solve. It Also happens to be the most wdl known eombinatoidAl optimization pacoblem, xad 
has Attracted interest ttcm n wide range of fidds. In operations research, mathematics and 
computer science, researchers have concentrated on algorithmic aspects. A particnlar fbcns has 
been on henristic algorithms — algorithms which do not gnarantoe optimal tonrs — for cases 
where exact methods are too slow to be of nse. The most effective henristics are based on local 
search methods, which start with a non-optimal tonr and iterativdy improve the tour within 
a well-defined 'iidghborhood"; a famons euunple is the lin-Kemighan henristic [l]. More 
recent eOorts have involved combining local search and non-deterministic methods, in order to 
refine henristics to the point where they give good enoi^h solntions for practical purposes; a 
powerful such technique is Chained Local Optimization [2]. 

Over the last fifteen years, physicists have increasingly been drawn to the TSP as well, and 
particularly to atockoAtic versions of the problem, where instances are randomly chosen finom 
an ensemlie. The motimtion has often been to find properties applicable to a large class 
of disordered systems, either throngh good approidmatc methods or through exjict analytical 
approaches. In our work, we consider two such stochastic TSPs. The first, the Euclidean 
TSP, is the more dassic form of the problem: N dties are placed randomly and independently 
in a ti-dimensional hypercube, and the distances between dties are defined by the Euclidean 
metuc. The second, the random Unk TSP, is a related problem devdoped within the context 
of disordered systems: rather than spedfying the positions of cities, we !g>ccify the lengths lij 
separating cities i and ^, where the Uj are taken to be independent, identically distributed ran- 
dom variables. The appeal of the random link problem is, on the one hand, that an analytical 
approach exists toi solving it [3,4], and on the other hand, that when certain correlations are 
neglected this TSP can be made to resemble the Endidean TSP. We therefore consider the 
random Unk problem as a random link appniximatian to the (random point) Euclidean prob- 
lem. Researchers outside of physics remain largely unaware of the analytical progress made 
on the random link TSP; one of our hopes is to demonstrate how these results are of direct 
interest in problems where the aim is to find the optimum Endidean TSP tour length. 

Our approach in this paper is then to examine both the Endidean problem and the random 
link problem — the latter for its own theoretical interest as well as for a better understanding 
of the Euclidean case. We begin by consideriqg in depth the Endidean TSP, including a review 
of previous work. We find that, given periodic boundary conditions (toroidal geometry), the 
Euclidean optimum tonr length Le averaged over the ensemlie of all possible instances has 
the finite size scaling behavior 



1 + 



Ficrn simulations, we extract very pacedse numerical lalues tat (iE(,d) at <f = 2 and d = 3; 
methodological and numerical procedures are detaUed in the appendices. We also give numer- 
ical evidence that the probability distribution of becomes Gaussian in the large N limit. 
In addition to these TSP results, we find a surprising universality in the scaling of the mean 
distance between ^th nearest ndghbors, for points randomly distributed in the tf-dimensional 
hypercube. Pinally, we discuss the expected behavior of /3e('^) ^ 'h^ large d Umit. 
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¥ig. 1. SeJf-flWHrapii^ af 3-D Eudidfian TSP aptdimim; oanvergienioe af ijfi(JV,3)/JV^-'^ an asequanioe 
□f raudaia iastanoea at increBaLii^ N. 



In Qic afcond part of the p^pcr ^ diaciBfi the r;mdiO(iii liiLt problem, confudcrlng it ;m 
Apprcadiiiation to the Enjdidican probkm. MAkin^ mm ot the cavity method, we mmpare the 
rADidom Imt /3rl(<^ with the Enclidejm iic (d) vaIuca obtjuncd team am siimiilAtioiB. Vic find 
that the rjuidom liiLt approadm^ition is eoxreet to within 2% At d = 2 .md 3. The reat of the 
flection iitndiefl the large d limit of the rjuidom Unt model ;ind its- implicAtionii for the EncUdean 
TSP. We aianunc analytically how ii^Lid) acalea at large and we relate the Ifd coefficient of 
the afiaodatcd poner aeriea to an underlying if-independent "renormaliBcd" modeL Finally, we 
present a theoretical analyfiia based on the Lin-KemigLan henriatic, anggrating strongly that 
the relative difference bet^en ,^)kl(^ and is positive andof f?(l/d^). The random Unt 

resnlts then lead to onr large d EncMean conjecture: 



(2) 



where 7 is Enler's constant. 
Z. Tha Euididiean TSP 

2.1. SCAIIKG AT Laejge N. — One of the caidicst analytical resnlts for the Endidean TSP 
is dne to Bcardwood, Halton and Hammerslcy [5] (BHH). The anthors considered N cities, 
distributed randomly and independently in a d-dimensional volume with distances between 
dtics given by the Euclidean metric. They showed that, when the volume is the unit hypercube 
and the distribution of cities uniform, Xe/JV^~^^^ is self-averaging. This means that with 
probability 1 , 

where (izld) is independent of the randomly chosen instances. This property is illustrated 
in Figure 1. (In &ct, the BHH result is more general than this and concerns an arbitrary 
volume and arbitrary form of the density of cities.) Eor a physics audience this large N limit 
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k cqnimknt, in AppsopuAtc nnita, to an infiiutc vdiuiciic limit At confitAiit dciuity. £e/JV^~^^^ 
thm concapoiiida to energy density th;it in sidf-Avciaging jmd i^a a well-defined infinite 
volnme limit. The original proof by BHH in quite complicated; simpler proo&i hxvc aince been 
given by Karp and Steele [6,7]. 

One of onr goak iA to determine (is (d). BHH gave rigoronfi lower and npper bonndfi an a 
function of dinujuion. Ear ai^ given iiuitancc, a trivial loner bound on Xe ifi the anm over all 
dtieai of the distance between i andita nearest neighbor in npace. In fact, siince a tonr at beat 
Unkfi a dty with itfi ifuo neaieat neighbors, thia bonnd can be improved npcn Iqt anmmii^, 
over all i, the mean ati'n nearest and nexLt-neareat neighbor diatancefi. Tiddng the eniunuble 
average of thia qnantity (that isi, the average over all inatancca) leadfi to the beat analytical 
lower bound to date. Ehr upper boundfi, BHH introduced a henriatic algorithm, now known 
afi "atrip", in order to generate near-optimal toura (diacusfied ako in a paper liy Armour and 
Wheeler [&] ). In two dimenaioufi the method involves dividing the square into at^acent cdnmnfi 
or strips, and sequentially visiting the dties on a given strip accordiqg to thdr positions along 
it. The respective loner and iq>per bounds give < (2) < 0.9204. 

In addition to bounds, it is possible to obtain numerical eattma-ttn fer BHH used 

two instances, N = 202 and N = 400, ftam which they estimated /)e(2) ^ 0.T49 usiqg hand- 
drawn tours. Surprisingly little has been done to improve upon this value in two dimensions, 
and essentially nothing in higher dimensions. Stdn [9] has found /3e(2) as 0.7G5, which is 
finequently cited. Only recently have better values been obtained, but as they come from near- 
optimal tours found by heuristic algorithms, they should be considered more as upper bounds 
than as estimates. Using a local search heuristic known as ''S-opt" [lO], Ong and Hnang [ll] 
have found ,^(2) < 0.743; using another heuristic, ''tabu" search, Elechter [12] has found 
/)e(2) < 0.731; and using a variant of simulated annealing, Lee and Choi [13] have found 
/)e(2) < 0.721. In what follows we shAU show what is needed for a more predse estimate of 
iizld) with, furthermore, a way to qnantify the associated error. 

2.2. EXTEACTIKG liEid). — As JV -v 30, L-e/N^'^/^ converges with probability 1 to the 
iustance-independent ii^id). Our estimate of (i^ld) must rest on some assumptions, though, 
since only finite mines of N are accessible numerically. Note first that at values of JV where 
computation times are reasonable, Xe has snhitantial instance-to-instance flnctnations. To 
reduce and at the same time qnantify these flnctnations, ne average over a large number of 
instances. ^Nc thus consider the numerical mean of Xe over the instances sampled, which itself 
satisfies the af^mptotic relation (3) but with a smoother convergence. Tb extract /3e(^)) 
must understand predsdy what this convergence in N is. 

If dties nere randomly distributed in the Iigrpercnbe with open boundary conditions, the 
dties near the boundaries would have fener ndghbors and therefore lengthen the tour. In 
standard statistical mechanical systems at constant density, boundary effects lead to corrections 
of the form surface over volume. Ebr the TSP at constant density, the volnme grows as N and 
the snrfiice as AT^'^^''. In a ti-dimeusional unit hypercnbe, then, the ensemble average of Xe 
would presumably have the large N behavior 

In order to extract ii^id) numerically, it would be necessary to perform a fit which includes 
these corrections. A reliable numerical fit, honever, must have few adjustable parameters, and 
the slow convergence of this series would prevent us £tom extracting (i^ (d) to high accuracy. We 
therefore have chosen to eliminate these boundary (surface) effects I^ nsii^ periodic boundary 
conditions in all directions. This should not change ,^(ii), but leaves us with fower adjustable 
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pjuramctcH ;md a. fnAtat convcxgciux, cnjiliiDg ha to watk. with auciJillcr valiua of N whicic 
numerical suniiiLitionfi xtc not too silow. 



Ebr th£ Ii^pcinnibc with periodic hoiLiiid;iry conditioiH, let na iabtodiiice the notation 



jVi-i/j ■ 



(5) 



where {Lj^) is. the average of £e over the eniunuble of infitancefu (/^^(JV,^) ia, in pliyfikal nnitji, 
the zero- temperature energy density.) We then wliih to nnderjitand how (is (N^ d) converges to 
its laige N llmit^ ,^('^)' standard statistical mechanical systems, there is a characteristic 
correlation lei^th Away &om a critical point, ^ is finite, and finite size corrections decrease as 
Q-^fi^ where W is a measure of the system "Vidth" . At a critical point , ^ is infini tc, and fini te 
size corrections decrease as a power of 1/W. Ebr diAenrdend statistical systems, however, this 
picture must be modified. Even If ^ is finite for each instance in the ensemble, the flnctnating 
disorder can stiU give rise to power-law corrections Sai ensemble averaged quantities. In the 
case of the TSP, this is particularly clear: the disorder In the positions of the cities induces 
Large finite size effects even on simple geometric quantities. 

To see how this might affect the convergence of liE{Nfd)f consider the following. Bex a 
given configuration of N points, call Di,{Nfd) the distance between a point and its Arth nearest 
neighbor, where k = 1,...,JV — 1. TaliL the points to be distributed randomly and nnifbarmly 
in the unit hypercnbe. Let us find {Ds,{Nfd)). Under periodic boundary conditions, the 
probability density p{I) of finding a point at distance J £Dom another point is simply eqiial (for 
<I < 1/2) to the surEice area at radius I of the <i-dimcnslonal sphere: 



(6) 



The probatiUty of finding a point's Jtth nearest neighbor at distance I (see Fig. 2) is equal to 
the probability of finding ^ — I (out of — 1 ) points within one point at I and the remaining 
N —k — 1 points beyond I: 



-k)pil)\l-j'f,il')di' 



N-b-i 



giving the ensemble average 

iD^iN.di} : 



r 1 



r(d/2+i) 



N-b-l 



di + 



(9) 



where the corrections are due to the I > 1/2 case, and are eifiponentiaUy small in N. 

Recognizing the integral, up to a simple change of variable, as a Beta function (Bla^h) = 
jji ^-i ^1 _ t)*-^ = T(a)£(h)/r(a -\- b)) pins a further remainder term eifiponentlally small in 
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k - 1 j>jinlB 



N - k- \. j5oinlB 



Fi^. 2. A pcnnt's ^ — 1 med^Hbaia: — 1 neaneat mei^lLbara are withia distance I, feth neareat 
njBL^bar ia si: 1, ami Temainui^ N — h — 1 pciLiita are beyand. i. 



JVj WC 9CC tluit 



r(J/2 + l/^^ r(ji:+ 1/d) T{N) 

^ r(Jt) r(jv + i/d) 



r(i) 



2Ar 



(10) 



Wc mnc con£noiitcd here with a mujirli^hlc, and hitherto micxplorcd, nnivcrflAlity: the tx&cst 
suuiQiC 1/JV iiciic!i giwa the JV-dcpoidiCiiee rcgardkaaof Jb. The sumie finite aise sifAliiig bohavior 
thiCiefjic Applira to all Jtth ncAieat ndghbox difitAiu:eA. 

It mi^t be hoped then th;it the ty]ic;il lint length in optiumni tom» wonld luive thifi N- 
dependence^ and that ^[N^d) wonld therefore have the aame 1/JV expanfiion. This ifi not 
qnite the case. The lint between dtiea i and j fignrea in the average (Da(JV, d)) whenever j 

the Jtth neighbor of it fignrea in hoirever, only when it bdonga to the optimal 

tour. Two diOcrent tinda of averages are being taten, and ao finite siize corrections need not be 
identical. Nevertheless, it remains ]iansible that ^[N^d) has a 1/JV series expansion, albeit 
a difierent one from (11 While we cannot prove this property, it is confirmed by an analysis 
of onr nmucrical data. 

Onr approach to finding ^[d) is thns as £}Bows: (i) we consider the ensemble average {Xe), 
rather than £c for a given instance, in order to have a quantity with a well-defined dependence 
on N\ (ii) we nse periodic bonndary conditions to eliminate snr&ce effects; (ili) we sample the 
ensemble using nnmcrical simulations, and measnre ,&(iV, d) within well controlied crrois; (iv) 
we cxiract ^{d) by fitting these values to a YjN series. 

2.3. Finite Size Scaling Heetlte. — Let us consider the d = 2 case in detail. found 
the most effective numerical optimization methods for our purposes to be the local search 
heuristics Lin-Kemighan (LK) [l] and Chained Local Optimization (CLO) [2] mentioned in 
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Fig. A. Finite aixe defpendeiiioe af tiio reacalad EudidiBaiL TSP aptimum. Best: S.t (x^ — Ji.Jiti) 
gi™: ,3jc(JV,a)/[H- l/(aJV)+ ■ ■ ■] = (\.nii(\{l-(\.(\lfl/N - 1.04a/JV^). Error bara repiBaant atatiatical 



the intiodiiiCtioii. Both hciudfiticsi, by dcfimtioa, give tom Icngtlui UluI axc not Al^ysi optimal. 
Howcvcx, it iii not ncccafuiry tlut the optimmu he fbnnd 100% of the time: theie is- ^ibccPidy 
A !iignifie;mt atAtifitical enor Aiifiing ftam infitjmee-to-inakince flm:tn;itionfi, jmd ao a fiirther 
jiyiitematu: enor dne to non-oiptimAl toiufi is. AceeptAble 2ls long asi thia error is kept negligible 
eompaied to the atAtiiitieAl eixor. Owt noethodii, along with relevant nnmeidcal detail, are 
dificnfHed in the appendices. For the preacnt pnrpoiiea, let un nmply noention the general 
natnre of the two henrifiticfi naed. LK wor]ia by perbiming a "triable-depth" loeal aeardi, as 
dificniHed further in Section CLO works liy an iterative process comlining LK optimizations 
with random perturbations to the tour, in order to explore many different local neighborhoods. 
We nsed LK fai ''small" N valnea (N < I7)j averaging over 250^000 instances at each valne of 
JV, and we nsed CLO far large" N valnes (N =30 and N = 100), averaging over 10,000 and 
G,000 instances respectively. 

We fitted onr resnlting ii^lN^d) estimates to a truncated 1/N series: the fits are good, 
and are stable with respect to the nse of snb^-samples of the data. Bex a fit of the form 
(ki^Ai = /JE('i)(l + A/N + B(N\ we find /Je(2) = 0.7120 ± 0.0002, with )^ = 5.57 for 
& data points and 3 fit parameters (5 degrees of beedom). Onr error estimate fior iij:(2) is 
obtained 1^ the standard method of perfiorming fits nsing a range of fixed valnes for this 
parameter: the error bar id) .0002 is determined by the valnes of /)e (2) which make eiuxcd 
its original resnlt by exactly 1, i.e., making = ^"^^ i^^ 'his case. 

It is possible to extract another /)e(JV, d) estimate by making direct nse of the nnivcrsality 
discussed ptevionsly: the nniverfial IfN series in (ll) suggests that there will be a faster 
convergence if we nse the rescaled data /JeC-W>2)/[1 -|- 1 /(SJV) + - This also has the appealing 
property of leading to a function monotonic in JV, aa shown in Figure 3. We find 



with the leading term having the same error bar of ±0.0002 as befere. Note that the 1/N term 
in the fit is small — 2 orders of magnitude smaller than the leading order coefficient — and so 
to first order the 1 + 1/SJV + " ' series is itself a good appacoodmation. 



errara. 




(12) 
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Fi^. 4. Difltributian of S-D Euclidfian TSP acaling variable Xjt = (ijc - {Lis:'j)/N^^~^^ . Shaied 
i^Span ia far N = H (llOOjQOQ ioataiiiCBa uaed) and aolid liae ia far H ~ AO (lOflOO instaiLoea uaed). 
Superimpaaed curve abawa (estrapalated ) Uioitiii^ CauaaioQ. 



The flame nicthiodology "stpa applied to the d = 5 case. The x^'n Again confLrmcd the fnnc- 
tiioaal tcntm cf the fit, and ^ find bom oni data ,^(3) = ± 0.0002. Alao, since csai 

initial work [l4], Johnfum e£ oJL have pcrfbiiiaed siiiiinlationa At d = 2, 3, 4, obtaining reanltn [15] 
eonaiatent with anin: /li:(2) ^ 0.7124, /)e(3) ^ 0.G9SD and /Ie(4) ^ 0.7234. 

2 A. Detribitidn of Optimum Tour Length g. — While BHH aM others [6, 7] have 
fihown that the variance cf Le/JV^~^^*' goea to zeio aa JV — v oo (aec alao Fig. 1), they have 
not detemoined how faat thia variance decicaaea. More generally, one might aali: how the 
distfibfitian at Lj:/N^~^f'^ hchAvca aa JV — i' so. Vfc axe awAie of only one leanlt, Ysy Rhcc and 
Thlagrand [IG], showing that the probaliUty of finding with \Lj^— {L^)] > t iaamaUer than 
KcTqt(—'^fK) for aome K. Unbrtnnately thia ia not atrong enough to give bonnda on the 
variance. 

Let na charajctciize the diatiibution at d = 2 by nnmerical aimnlation. Ebr motivation, con- 
aider the analogy between Le/JV^~^^^ and E/V^ the energy denaity in a diaordeied atatiatical 
ayatem. If the ayatem'a correlation length ^ ia finite (the j^atem ia not critical), E/V haa a 
diatribntion which boconica Ganaaian when V — ao. Thia ia becanae aa the anbvdnmea in- 
creaae, the energy denaitiea in each anbvolnme become nnicorrelated; the central Umit theorem 
then appUea. A conaeqnence ia that ^r^, the variance ofE/V^ decreaaea aa . If ^ ia infinite 
(the ayatem ia critical), then in general the diatribntion of£^/Vianot Ganaaiam In both caaea 
though, the aelf-avcraging otE/V anggeata that the acaling variable X = (E— {E})/<tV haa a 
limiting diatribntion when V — so. 

In the caae of the TSP, it can be argned naing a theoretical analyaia of the LK henriatic that 
at > 2 the ayatem ia not critical. By analogy with E/V^ if we talie anbi^olnmea to contain a 
fixed number of citica, the central limit theorem then auggcata that L-E/N^~^f^ haa a Gauaaian 
diatribntion with decreaaing aa JV"-'. The acaling variable Xfi — (it — (iE})/JV^^^~^^^ 
ahould conaequently have a Gauaaian diatribntion with a finite width for AT —t- so (and at 
<f > 2). Numerical reaulta at d = 2 (aee Pig. 4) give good anpport for thia. 
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CoNJEcrruHES OK TEE Laeuc^e d Limit. — In mcs^t sitAtkticAl ucuoclLaiiksi paDolicmfi, the 
krgc dinoiCiHioiuiL limit iatcodiuxa almplificatlonfi bccanac fliiftiuitionfi hccumc ntgli^blcL Bat 
the TSP, can one cxipcct (iEi^i to Lave a jiimple limit aa d — i- so? AgaLn, confdder the pBopcxty 
of the kih neareat neLghbor distance Db. In the large N limit, (ll) givesi 



NiDi,{N,d)) N' T{k) > ™ ^' ^'^g^ 



(13) 



1 + ^+ (14) 



where = — y -|- + + - - - (7 ifl Enler's eonatant). Notice that ^ Ini: at large Thia 
anggeata atrongly that nnleaa the "typical" Jtnsed in the optimum tonr growa eifiponentiaUy in 
we may write fcr d — ^ so: 

Up to Q{\ld)^ thia e]q>rea8ion ia identical to the BHH lower bomid on ,^(^) diacnaaed in 
Section 2.1, given by the large limit of JV^^^(£ii(JV,d) + £»s(JV,d)}/2. 
A wealier coiqectnre than (15) haa been patopoaed by Bertaimaa and van Ryzin [l 1^: 

tit:{d) -J ■/d/2iw aa d -^ 30. (16) 

Thia limiting behavior waa motivated by an analogona leanlt br a related combinatorial op^ 
timizatlDn problem, the minimnm apanniog tree. Unfortnnately, there ia no proof of either 
(15) or (IG); in particular, the upper bound on ,^(^) given by atrip, diacuasied in Section 2.1, 
behavea aa y/d/G at large d. Thua if the conjecturea are true, the atrip conatmction leada 
aaymptoticaUy to toura which are on average 1 .69 timea too long. Can we derive atronger 
upper bounda? A number of heuriatic conatructlDU methoda ahould do better than atrip, but 
there are no reliable calcnlationa to thia effect. The only improvementa over the BHH reaulta 
are due to Smith [IS], who generalized the atrip algorithm 1^ optimizing the ahape of the atripa, 
leadiug to an upper bound which ia \/2 timea greater than the predictiona of (15) and (IG) at 
large d. 

In apite of our inability to derive an upper bound which, together with the BHH lower bound, 
would confirm the two conjecturea for <f — so, we are confident that (15) and (IG) are true 
becauae of non-rigoroua yet convincing argnmenta. One ia a proof that (IG) ia aatiafied for the 
TSP if it ia aatiafied far another related comlinatorial optimization problem (aee Appendix: 
D tai detaila). A more powerful argument, pieaented in Section 3.G, reUea on a theoretical 
analyaia of the LK heuriatic. It auggeata that up to 0{l/iP)f fh:id} ia given by a random Hnt 
approximation, leading to a conjecture even atronger than (15). 



3. The HBndoiu Link TSP 

3,1. CoMLEgpoKEEKCE WITH TEE EccuDEAN TSP. — Let ua iKw conalder a problem at 
firat aight dramatically different bom the EucUdean TSP. Inatead of taMug the poaitiona of 
the N dtiea to be independent random variablea, talu: the leugtha I^j = Iji between citiea i 
i (1 £ £ -^) to be independent random variablea, identically diatributod according to 
jiome We apeal;: of leqgtha rather than diatancea, aa there ia no diatance metric here. Thia 
problem, introduced hy pliyaiciata in the 19S0a [19, 20] in aearch of an analytically tractalie 
form of the traveling aaleaman problem, ia called the vnnd&m iink TSP. 
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Th£ coniifftion hctwcf n thk TSP and the Em:lidcAn TSP k not obi^ioiH, wc now hxvc 
rADidiOm liuka lathcr tLau laiiidDiii pcdntn. NcvKxthidcfifi, one c;m relate thf two probkuoa. To 
fifc thia, cuiHidcr the probability diatiibntion for the difitjuiice I between a fiiu:d pair of dtiea 
(ij) in the Enclidean TSP. Thifi diatiibntion, in the nnit hypermbe with penodic boundary 
conditionfi, m given for < I < 1 /2 by the cxpieaflion in (G): 

Of conrfie, in the Endidean TSP the lint lengthfi are Yiy no meana independent random vari- 
ables: cunelationfi hmcJi aa the triangle ioeqaality are present. Honever, aa noted by Ml^zard 
and Pariai [3], cunelationfi appear eud naively when eonaideiing three or more diatancea, ainee 
any two Ehiflidean diatanicesi are neceafiarily independent. Let mi adopt as. the diiitribn- 
tion in the limit of amall I for the random link TSP, where d in thia case no longer represicnts 
ph^aical dimension bnt is simply a parameter of the model. The Enelideaii and random UnJi: 
problems then have the same small I one- and two-link distribntions. In the Large N limit the 
random link TSP may therefboDe be considered, rather than as a scpATAtc problem, as a random 
link approxifnation to the Endidean TSP. Only joint distubntions of three or more links differ 
between these two TSPs. If indeed the corrdations involved are not too impartant, then the 
r;indom link tan be taken as a good estimate of (i^ld). We shall see that this is tme, 

particnlarly for large d. 

3.2. SCAUKG AT Larc^E N. — As in the EncMean case, we are interested in nnderstanding 
the N iX) scaUi^ law in the random link TSP. It is rdativdy simple to see, following an 
argnment similar to the one in Section 2.2, that the nearest neighbor distances have a 
probability distribation with a scaling £\ctor JV"^^** at large JV. Vannimenns and Mdzard [20] 
have suggested that the random link optimnm tour lepgth with N links will then scale as 
JV^"^^'', and the tonr will be sdtaveraging, i.c. 



parallel to the BHH theorem (5) for the EncUdean case. This involves the implicit assumption 
that optimnm tours sample a representative part of the Dt, distribution, so no further N 
scaling effects are introduced. The assumption seems reasonable based on the analogy with 
the Euclidean TSP, and font our purposes we shall accept here that /3rl(^) exists. However, 
there is to our knowledge no mathematical proof of self-averaging in the random link TSP. 

EloUowing the discussion of Section 2.1, let us consider some bounds on the enscmlic average 
{^^rl} ^ derived in [20] . As before, we get a lower bound on (iRLid) using nearest and next 
nearest ndghbor distances. Eor an upper bound, the "strip" algorithm used in the Euclidean 
case (Sect. 2.1) cannot be applied to the random link case. On the other hand, Vannimenns 
and M^zard make use of an algorithm called "greedy" [21] : this constructs a non-optimal tonr 
by starting at an arbitrary dty, and then successively picking the link to the nearest availalie 
dty until all cities are used once and a closed tour is formed. At d > 1 , greedy gives rise to 
tour lengths that are sdf-averaging, and leads to the upper bound [20] 

i/ir a — 1 

At d = 1, the presumed scaling (IS) suggests that is independent of JV, whereas greedy 

generates tour lei^ths which grow as InJV. There is numerical evidence [4,22], however, that 
the if = 1 model does indeed satis^ (l&), and that ,^(l) w 1.020S. 
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3,3. SoLcnoK wia TEE CAVITY EquAHOKg. — Slmx the wor];: of VAimiiiacaiifi And M^aid, 
fifvcxAl giDoapA [23-25] hxvc tucd to ''jioLvc" the atAtkticAL iiQ£diADic;il piohkm of the r;iiiidom 
link TSP At finite tcmpcrAtnuc vaing the icplicA method, a tedmiqnc developed for jmaljrzin^ 
difioideitied syAtcnaa. smdh sipin ghiasiLS. [26]. Tb date, it h^ji only been poafiible to obtain part 
of the high temperature aexiea of thisi siysitem [23] . In view of the intractability of thcfie replica 
approaches, Mfzaxdand Parifii have derived an analytical fioilntiDnnfiing Another technique &om 
spin gLafm theory, the "cavity method". The details of thia approach are beyond the acope of 
thifi paper, and are dificnfiaed in several technical axtidea [3,2G,27] . Eat readers acquainted with 
the langnage of difiordered aysitemfi, hoirever, the htoad outline ia afi blloiriii: one beginii with 
a representation of the TSP in tcrma of a Heificnbcrg (mnlti-dimensional apin) model in the 
limit where the spin dLmenfiion goes to zero. Under the assumption that this system has only 
one cqmlibrium state (no replica symmetry brealdng), Mfsaid and Paiisi have then written a 
recursion equation fbor the system when a new (N + l)th spin ia. added The cavity method 
then supposes that thia new spin's cOcct on the N other spins is negligible in the Lai^e N limit, 
and that its magnetization may be expressed in terms of the magnetisations of the other spins. 

Using this method, Krauth and Mt^zaid have derived a self-consistent equation fer the ran- 
dom link TSP, at — i' [4] . They have determined the patobabiUty distribution of lint lepgths 
in the optimum tour in terms of ^ii(jr), where is the solution to the integral equation 



These equations can be solved numerically, as well as analytically in terms of a 1/d power series 
(see next section). At d = 1, Krauth and 3Ji(!zard compared their parediction with the results 
of a direct simulation of the random Unt model; their numerical study [4,22] strongly suggests 
that the cavity prediction is exact in this case. It has been argued, furthermore, that the 
cavity method is exact at — i- for a-nj} distribution of the independent random Unks [2&] . 
Good numerical evidence has been found fcx this, notably in the case of the matching patoblem, 
a related combinatorial optimization problem [26i] . The validity of the cavity assumptions 
there&ie does not appear to be sensitive to the dimension and ^ shall assume that (21) 
holds for the random lint TSP at all d. 

Krauth and Ml^zard computed the d = 1 and d = 2 cases to give ,^l(1) = 1^(I2DS and 
/1kl(2) = 0.T251. Since (inLid) is tatcn to approximate /1e(^, let us compare these values 
with their Euclidean counterparts. At d = 1, the EucHdean TSP with periodic boundary 
conditions is trivial (,^(l) = 1); the random lint TSP thns has a 2.1^ relative excess. At 
d = 2, comparii^ with ,^(2) = 0.7120 found in Section 2.3, the random lint TSP has a 
excess. In low dimensions, the random lint results aie then a good approximatiou of 
the Euclidean results. The approximation is better than Krauth and Mf^zard believed, since 
they made the comparison at d = 2 nsiug the considerably oi^erestimated EucHdean lalue of 
/Je(2) as 0.749 &om [S]. 

Extending the numerical solutions to higher dimensions, at d = 3 we find /^hl(3) = G.71D0, 
which compared with /Se(^) = 0.6979, has an excess of 1.7^. Some further random lint values 
aic /haX'i) = 0.7322 and ,;Jkl(5) = 0.7639. The value at d = 4 may be compared with the 
EucHdean estimate of Johnson est ai [15], ,^(4) ^ 0.7234, giving an excess of 1.2^. The /Se(^ 
data at d = 1,2,3,4 Ihcrcfcrc si^gest that the random lint approximatiou improves with 
increasing dimension. This leads us to study the Hmit ^dien d becomes large. 





Their probability distribution leads to the prediction 




(21) 
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Fig', ii. riiiTMiTuiifiTiAl dependence af Teacaled Tandom. link TSP aptuDoim, ahiown by amall paints, 
betTween Danwrgin^ "^eedy" upper bound (dotted line) and neaneat-nei^hboEa Iowt bound (dashed 
line). Plua ai^a at ^ = :S and d = akow Eudidiean reaulta Iot comparidan. 



3.4. De^ensional Dependence. — The Ijugc d limit w^si coiuidcicd VjumiiiQciuifi juul 
Wz;ird [20]. Ebr /3m,(d), Ihc lower hound obtjuncd feom {B^iN^d) ■^D^{N^d))ll. by way of 
(11) and the iq>pcr boimd given in (19) diBcr at large d only by 0(l/d), giving: 



l+o(l)]. (22) 

Note that thia enaci reanlt i» the random lint analogue of the Euclidean conjectnre (15). 

Ebr valnca of d < 50, wc have caknlatcd ^mXd) numerically naiug the cavity equations 
(20,21). The rcfinlta are ahiown in Elgnre 5, along with the convergiDg iq>per and lower bonnds, 
and onr low d Euclidean reaultsi. 

Eor large we may siee whether the cavity cqnationa are compatible with (22) liy aclving 
them analytically in termsiofa 1 power aeriea. Define Qi{x) = Qd{C(d+ l)^'*" [l/2-|- x/d]). 
(20) may then be written: 



,ftii,(ff) 



= — 



(^d) 



dy 



(23) 

-Oiiv-i^y. (24) 



Strictly apeaJdng, the expansion of (l + [i -|- y]/dy^ ^ ia only vaM in the interval —i—d< 
y < —x + d; however, for large y it can be shown that Qi{y) ^ j/^, ao the e"^*"'^^ term in the 
integrand maJiea the j/ > —x + d contribution exiponeutLaUy amaU in d. 

Furthermore, extending the integral 'a lower limit to iuclnde the region y < —x — d alac 
contribntea a remainder term enponentiaUy amaU in d. If we write the integral with ita lower 
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liimt at = —SO) the cqnatioa may he solved: 

[ i V 2 2 S J 

+°{^)]' 

whiCic 7, wc recall, rcpccsunitA Enkr'a conatant. Usiiivg (21), wc thm find 

whiidi k perfectly compatible with (22). Thifi provider fiutther evidence that the cavity uiethod 
ifi e^Lact for the random link TSP. 

3.b. Renoemauzed ILAKEOM Likk Model at Laese d. — We can motiiate the large 
d acaUqg found in the panevionfi section by examining a diflerent sort of random link TSP. 
Consider a new '^noimalizcd'" model where Unt lengths'" i^j are obtained &om the original 
lij by the Unear transformation = d[Iij — {Di{N^d))]/{Di(N^d)). Note that the jr^ may 
talie on negative valnes, and that the nearest neighbor length in this new model has mean zero. 
Since the transfbimation is Unear, there is a direct eqnivalence between the renormalizod x^j 
and original Ijj TSPs, and the two have the same optimnm tonrs. The renormalized optimum 
tonr length Lj, may then be given in terms of the original tonr length Li Yiy 

L,-NiD,iN,di} 
^'-^ {D,{N,d)) ' f^^J 

Now taJie N and d sa. It may be seen finom the Uj distribntion (l7) and the 

{Di.{N^d)) expansion (14) that the random varialies Hj have the d-independent probability 
distribntion p{x) ^ e3tp(x — 7). Also, in the large N limit, since Lf scalesas JV^"^^*' and 
{D^) scales as N~^^^^ we c]:pect {LJ) Nji for some ji which mnst be, like /)(x), independent 
of Then, from (27), the TSP in the original lij variables satisfies 



{^,) ^ N{D,{N,d)) (2S) 
or, using the e^qyansion (14), 

This lesnlt may be compared with onr cavity solution of (2G), where the \jd coefficient is eqnal 
to 2 — la2 —27. If the cavity method is correct at 0{\jd)^ which we strongly beUcve is the 
case, then a direct sdutiDn of the renoxmaHzed model should give = 2 — lu2 — 7. Wort is 
cnrrently in pcogress to test this claim by numerical methods. 

3,G. Large d Acclhacy of tee Random Link Approximation, — Since the random 
Unt model is considered to be an approcdmation to the Euclidean case, it is natural to ask 
whether the approximation becomes exact as d — v ao. In this section we argue that: (i) in 
stochastic TSPs, good tomrs can be obtained using almost exclusively low order neighbors; 
(ii) the geometry inherent in the EhicHdean TSP leads to l^{d) < f^iaJid) in all dimensions d\ 
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Fi^. B. RecuraivB ocaLstructdciiL of lewavei Uiika (daabed lines) and added Unka (bold Unea) ia an 
LK SBaich. 

(ill) th£ icLitivc cnor of the r;iiiidiOm liuik AppDoadniAtiDa docicjifica ^ At Large d. All tlucc 
djuiufi Aic hiAfifd on a thjcorctical AnAlyaifi of the Lia-KcmighAa (LK) hcmifitk algonthm fat 
mmtmctiii^ ncar-optiiiCLal toim. 

The LK Algonthm "soTtiA oa faHowa [1,39]. An LK Atanh at Arts with xd. aihitixty tour. The 
priiiiciplc of the search k to nnhiAtitnte Unka ia the tour recnraivdy, as- UlniitrAtcd a£hiemAtk;illy 
in Pignre G. The firat atep eonaiata of choofiiDg .m jirbitraiy fitjirtiiig dty ia. Call ii the next 
dty on the tonr, And fi the link between the two. Now remove thda link. Let he the nejueat 
ndghbor to ii th^t was not connected to ii on the oiiginAl tour, xnd let J' he a new link 
connecting ii to i\. We now h^ve not a tonr bnt a ''tAdpole graph", contmning a loop with a 
t;ul Attached to it at At thda point, caU one of the citiea next to i\ on the originAl tonr, 
And remove the link 1^ bet^en the two. There Aie two posHibilitiea for (and thna fg): LK 
choffiua the one which, if we wote to pnt in a new link between And to , wonld give a dngle 
doaed tonr. Now aa befeie, Let be the neAreat neighbor of that WAa not connected to 
on the originAl tonr. And let 1'^ be a new Unk bet^en the two. Thia gives a new tAdpoLe. The 
proceaa continnea recnraivdy in thia mAnner, with the vertex hoppii^ Aronnd while the end 
point s-tays. fixed, nntil no new tAdpdea Are fbnnd. At eAch atep, LK chooaea the new i^ao2A 
to aUow the pAth to be doaed np between And to , fbrmiitg a aingle tonr; the Deanlt of the 
LK aeATch is then the beat of aLL anch doaed np tonra. The LK al^riikm conaiata of repeAting 
theae LK aeArchea on different atArting pointa to, eAch time naing the cnrrent beat tonr Aa a 
atATting tonr, until no further tom improvemcnta Are poaaible. 

Let na fiiat aketch why the LK Algorithm leAda to tonra which nae only Unka between '^ncAr'" 
ndghbora, where 'iieAr" meAna thAt the ndghborhood order k ia aniAU And doea not grow 
with d. Conaider Any tonr where a aignificAnt &Action of the Unka connect diatAnt neighhora 
(Lutgc k). The Linka I'^ which the LK aeAich anbatitntea for the Are, Yiy definitiDn, bet^cn 
very neAr ndghbora (k < 3). Aa Loqg Aa mAuy Long Unka exiat, the probAliUty At eAch atep of 
anbatitntipg a neAr ndghbor in plAce of a &r neighbor ia aignificAnt. TbwAida the beginning of 
An LK aeArch thia probAbUity ia rdAtiveLy conatAnt, ao the expected tAdpole length will decreAae 
UncArly with the number of atcpa. Even tAking into Account the fact thAt closing up the pAth 
between t^ And t^ might requite inaerting a Unk with j: > 3, there ia a high probAbUity Aa 
N iXi thAt the improi/ement in tAdpole length fAr ontwdgha thia coat of cloaing the tonr. 
Thna tat atochAatic TSPa, regArdleaa of if, the LK Algorithm cau At lArge N replAce a11 bnt a 
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tiny fraction of the ko^ liuikfi with ^hoitt liiLl». It fbLlowa thj\t in. ju:aii»l;uu:c with awt Eiiiclidicaii 
TSP ASHmuptioQ at Scctioa 2.5, the "typical" k nad in the optimimi tour rciiQAiiia!iiu;ill At large 
d. Thifi providca very powerful finpport for the /S^ (d) cnnjcctnrica (15) and A conficqiieiux, 
ma]dqg mm of the cxLact juyiuptotic (hiiid) HAvlt (22), ia. tkat the relative dlBerenee bet^ea 
Ihiid) and ,iJni,(d) is at moat of 0(l/d). 

Onr second argument canoemfi why (iB.i,(d) mnat he greater than (izid) at all d. Bat the 
random lint TSP there ia no triangle inequality, whieh meam that given two edgea of a triangle, 
the third edge ia- on average loiter than it would be for the Eudidean TSP. Applying thia to 
our LK search, we can cxipect the Unt between and ia cloaing up the tour to be longer 
in the random lint caae than in the Euclidean caae. Thua on average, the LK algorithm will 
find longer random lint toura than Euclidean tonra. In fact, thia property hdda aa well foot 
any LK-Ute algorithm where the method of chooaing the and I'^ Hnta is. generalized. If the 
algorithm were to allow all poaaibiUtica toi I„i and 1^ , wc would be aurc of obtaining the eiuct 
optimum tour, given a long enough aearch. In that caae, the ineqnality on the tour lepgtha 
found by onr algorithm Icada directly to ,')hl(<^ > ,'^e(^)' aurpriaingly, the numerical 
data confirm thia inequality at d up to 4 (althoi^h one ahould be cautiona when applying the 
argument at d = 1). Note alao that the inequality in itaelf impUea conjecturea (15) and (IG) 
for the Eudidean model, aince it auppUea predaely the upper bound we need on (icld). 

Finally let ua exijiain why the relative difierence between ,f)KL(d) and ,^(^) ahould be of 
0{l/iP). Thia involvea quantifying the tour length impatovement difumaaed aboi^e. It ia dear 
tLat any non-optimal tour can be improved to the point where Unta arc moatly between 
ndghboota of low order. If LK, or a generalized LK-Ute algorithm, ia able to improve the tour 
further, the rdative diBerence in length will be of f}(l/d); aee thia from (14), noting tliat the 
neighborhood order k ia amall both before and after the LK aearch. Now we need to qnantify 
the prcb&bility tLat LK indeed aucceeda in improving the tour. We may conaider the vertex of 
the LK tadpole graph aa exjecuting a random wait, in ^lich case the probability of doaing up 
a tour liy a auffidently ahort Hut ia equivalent to the probability of the random wait 'a end- 
to-end diatance being anfficiently amall. In that cafie it may be ahown that, over the courae of 
an LK aearch, the probaUUty of ancceaaMly doaing a random Hnt tour minna the probability 
of aucceaafuUy cloaing a EncMean tour acalea at large d aa 2/{d — 2). Ftam thia, we conclude 
that improvementa in the Eudidean model are 0{l/d) more probable than in the random Hht 
modd. Now, the rdative tour length improvement fer the Euclidean TSP compared to the 
random lint TSP ia aim ply the rdative tonr length impacovement taken & better towr iafotind^ 
timea the probability of findiitg a better tour — hence 0{l/d^), If we conaider a generaUzod 
LK aearch aa deacribed in the previona paragraph, where the algorithm neceaajuily finda the 
tme optimum, then thia reault appHea to the exact (i'n: the relative difference betiseen li^Li^ 
and (izidi wUl acale at large d aa l/tP. 

Three commenta are in order concerning thia anrprifiingly good accuracy of the random Uht 
approximation. Eirat, the factor 2/(d — 2) ia only appropriate for large d. It ia not amall even 
for d = 4. (Ita divergence at d = 2 ia aaaociated with the £ict that a two-dimenaional random 
wait retnma to ita origin with probability 1.) ^Nc therefore expect the 1/d^ acaHng to become 
apparent only for d > 5, beyond the range of our numerical data. Second, we have aeen that 
the coeffident of the 1/d term in iiari'^) be obtained by the cavity method. Aaanming 
that thia method ia correct and that .^^(d) and .^(d) do indeed converge aa 1/d^, thia leada 
to a particularly atrong conjecture for the Eudidean TSP: 
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Thiid, thin type of LK AnALysiia cxa in £u:t he oittndcd to mAay othiot mmhiiLatoml optimizA- 
tloa probkuofi, audi ^ the ;ifiaigimi(Jit, imtdiiiig And bipartite matdiing pcolicma. In thcac 
we e]q>cet the rjuidam liiLt approodmatioa to give liae to a Q{\jS ) rel^itive error jmt ra 
in the TSR 

4. S nrninaiy aul CoDjClllsiaDS 

The first goal in onr work hzA been to inveatig;ite the finite size acAling of Xc, the optimnm En- 
didean traveling aalefinian tonr lengthy and to obtain precise estimates fcr ita large N beliaviorL 
Motivated by a remarkable nniversAlity in the Jtth nearest neighbor distribntion, we Lave found 
tfiat nnder periodic boundary conditions, the convergence of {Xe^/JV^^'^^^ to its limit fiE{d) is 
described \^ a series in 1/^^ This Las enabled us to extract /3e(2) ^ind ,^(3) nsiug numerical 
simulations at small values of JV, where errors are easy to controL Furthermore, thanks to a 
bias^free variance reduction method (see Appendix: B), these estimates are exitremely precise. 

Our second goal has been to ouumne the random Link TSP, where there are no correlations 
between Unk lengths. We have considered it as an approadmation to the Euclidean TSP, in 
order to understand better the dimensional scaling dt jhiid). Fbr small we have used the 
cavity method to obtain numerical values of the random link /^KL(d). Comparing these with 
our numerical values fcr ^^{d) shows tfiat the random Hnk approximation is remarkably good, 
accurate to within 2% at low dimension. Fbr large d, we have solved the cavity equations 
analytically to give jisa,{d) in terms of a 1/d series. We have then argued, using a theoretical 
analysis of iterative tonr improvement algorithms, that the relative diQerence between .^)hl(^ 
and Ae. (d) decreases as . This leads to our conjecture (3/0) on the large dbehavior of (<f), 
specifying both its asymptotic firm and its leading order correction. 

Let us conclude with some remaiuing open questions. Pbrst of all, while the cavity method 
most likely gives the ex^ct result fer the random Unk TSP, we would be interested in seeing this 
argued on a more fundamental physical level. Readers with a bacliground in disordered systems 
will recognize tLat the underlyiug assumption of a unique equilibrium state is false in many 
MP-complete proliems, and in particular in the spin-glass problem tfiat has inspired the cavity 
method. What makes the TSP different? Second of all, our renormalized random Unk model 
provides an alternate approach to finding the 1/d coefficient of the power series in (iKL(d)f 
and could prove a usefiil test of the cavity method's vaUdity. A solution to the renormaUzod 
model using heuristic methods appears within reach. Third of all, the 0(l/^) convergence 
of the random Unk approumation merits further study, £Dom both numerical and analytical 
perjq>cctives. NnmericaUy, EhicUdean simulations at <f > 5 could provide powerful support for 
the form of the convergence, and thus for our conjecture (30). Analytically, the qualitative 
arguments presented in Section 3.G, based on the LK algorithm, could perhaps be refined by 
a more qnantitative approach. Lastly, it is worth noting that the 0(lf^) convergence should 
apply cqnaUy weU to the di/itribiiiton of link lengths in the optimnm tour. The random Unk 
prediction for this distribution can be obtained from the cavity method [4]; an interesting test 
would then be to compare it with simulation results for the true Endidean distributiDn. 

Acknowledgiiieiits 

We thank E. Bogomolny, D.S. Johnson, M. M^zard, S.W. Otto and N. Sourlas for fruitful 
discnssiions. NJC acknowledges a grant from the Human Capital and Mobility Program of the 
Commission of the European Communities; JEblM acknowledges finandal siq>port from the 
F^nench Ministry of Higher Education and Research; OCM acknowledges siq>port from NATO 
travd grant CRG 920631 and item the Institnt Universitaire de France. 



44 



Reprint: Journal de Physique 



N°l RANDOM LINK APPROXIMATION FOR TIIE ElIOUEEAN TSP 1,11 



Appendbc A 

Overview of the Nmuerical Methodology 

hi the fblkiwliig, wc dififiififi the proccdmca nacd to obtain the ixw data ttcm which (izid) 
and the finite aixc acahng cucffificnta juc cxtrajctcd. Two mAjor problaiofi mnfit he acUvcd in 
order to get good eatluQJiteji of liElN^d). Fiitfit, ii^iN^d) in defined as. An enfiemble Averjigc 
{LE(Nfd))fN'-~'-f'^f but ia mcMnicd by a numeric;!! Average over a finite aamjie of iufitanceA. 
The inatanee-to-iufitance fluctuatioufi in Lj^ give riae to a fitatifitical error^ which decreafiea only 
as. the inverse square root of the sample size. Keeping the statistical error down to acceptable 
levels could require inordiuatc amounts of computing time. Vfc therefore find it useful to 
introduce a variance reduction trick: instead of measuring £e, we measure Le — AL* , where A 
is a free parameter and L* can be any quantity which is strongly correlated with X^. Details 
are given in Appendix: B. 

A second and more basic problem is that it is computationally costly to determine the 
optimal tour lengths for a large number of instances, precisely because the TSP is an NP- 
complete problem. The most sophisticated "branch and cut" algorithms can talje minutes on 
a workstation to solve a single instance of size N < 100 to optimality. However, we do not 
need to guarantee optimalLty: the statistical error in /Se (JV, already limits the quality of our 
estimate, and so an additional (systematic) enor in is admissible as long as it is negligible 
compared to the statistical error. We may thus use fast heuristics to measure Xf, rather than 
exact but slower algorithms. This is discussed further in Appendix: C. 



AppQDdbc B 

Statistical IVroTS and a Variance deduction THck 

Consider estimating {Lj^(Nfd)) at a given N Yiy sampling over ma ny instanc es. If we have 
M independent instances, the simplest estimator for {LE{Nfd)) is Xi:(JV,d), the numerical 
average over the M instances of the TmiTiinmum tour lengths. This estimator has an expected 
statistical error <r(M) = 'S'l^^/^/M, where tTL^ ^ instance-to-instance standard deviation 
of Xe. 

Now let us define Xii to be the sum, over aU cities, of i:th nearest neighbor distances. {Xi,) 
is its ensemble average; in terms of the notation used earUer in the texi, {Ls,) = N{Di, ). It has 
been noted by Sonrlas [30] that X^ is strongly correlated with Xi, X^ and Xg. He therefore 
suggested reducing the statistical error in (Xe) using the estimator 

Es = (Xisa}XE/Xiaa, (B.l) 

where X193 is the arithmetic mean of Xi, X9 and Xg. The ensemble average {Li^i} ciui be 
calculated analytically ttcm (11), and so the variance of Es comes from fluctnations in the 
ratio Xx/Xiga. If X^ were a constant fictor times Xua, this estimator would of course be 
per±ct, i.e., it would have zero variance. This is not the case, however, and furthermore the 
use of a ratio liases the Somdas estimatort its true mathematical expectation value differs 
from (Xe(J\^, d)} Yiy 0(l/JV). To improve upon this, we have introduced our own bias^free 
estimator [31]: 

£^_p = A{Xis J -|- Xe — AXij J (B.3) 

where Lis ^ 'be arithmetic mean of Xi and Xg, and A is a free parameter. Our estimator has 
a reduced variance because X^ and Xu are correlated. It is easy to show that the variance 
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of £ii_p iii miniiuiBcd a.t a niiiqnfj valnc of A, A* = C(LEfLi2)<^Li,f'^i.:a) 'S'licrc C(AfB) = 
{(■A — {A)){B — {B))'^f<Tji<Tg is th£ CDUcLitiDii cocffidcut of A .ind B, The vjm;iiicc then 
hhccomis ^^L^.^. = ""Lu [l — (£ej iia )] Empiidcallyj wc have fonnd thia vaiianjcc reduction 
pioccdnie to he qiiite eflective, aince ^1 -C^ aj O^S at if = 2 and ^1 - ss 0.31 at i = 3. 
The atatiatical error k thnfi reduced by about a factor of 3; thifi meana that for a given error, 
compntiug time in reduced li^ about a factor of 10. 

Appendbc C 

Cbntnol of Systematic Eftots 

Our procedure for esitimati^g Lj^ at a given instance involves running a good heuriatic m 
timea ftcm random atarta on that inatauce, and taking the beat tour length found in thoae 
m tuala. The eifipectcd aysitematic error can be bund bom the fitfquenciea with which each 
local optimum appeara in a large number of teat triala. (Thia large number muat be much 
greater than the actual number of triala uacd in production runa.) The meaaurement ia 
perfomiod on a auffidently large aample of inatancea, finom which we extract the average adze 
of the ayatematic error in {Lj^^N^ d)) aa a fimctiou of ui. We have feund that in practice, thia 
atxor ia dominated by thoae iubcquent inatancca where a anVoptimal tonr ia obtained with the 
higheat frequency. 

Aa N increaaca, the probability of not finding the true optimum incrcaaea rather faat; for a 
given hcuriatic, it ia thua neceaaary to increaae m with N in auch a way that the ayatematic 
error remaiua much amaller than the atatiatical euor. If the heuriatic ia not powerful enough, 
III wUl be too large for the computational reaourcea. Bat our purpoaea, we have found that 
the LinrKemighan heuriatic [l] ia powerful enough for the amaller valuea of (JV < 1 7). For 
W < N < 100, it waa more cffidcnt to awitch to GLained Local Optimization (CLO) [2,32], 
a more poweifnl heuriatic which can be thoi^ht of aa a generalization of aimulatcd annealing. 
(When the temperature parameter ia aet to zero ao that no up-lull movea are accepted, aa waa 
the caai: for our runa, CLO with embedded Lin-Kemighan ia called Iterated LinrKemighan'" 
[33,34].) With theae chodcea, uaing in two dimenaiona m = 10 fat N < 11 (LK), m = 5 for 
N = 30 and m = Wfcac N = 100 (CLO), we have kept ayatematic erroxa to under 10% of the 
atatiatical erroxa. 

Appendiic D 

BouDidiDg iisi^ using the Bipartite Matdung Problem 

Given two aeta of podnta Pi , . . . ,P jv and Qi , . . . , Q jv in d-dimenaional Eudidean apace j the 
bipartite matchii^ (BM) problem aakafor the minimum matching coat L^u between the P^'a 
and the Qi'a, with the conatraint that only linka of the fatsn P — Q are allowed. The coat of a 
matching ia equal to the aum of the diatancea between matched paira of pointa. When pointa 
Pj and Qi are choaen at random in a d-dimenaional unit hypercube, it ia natural to expect 
ZiHii/-W^^~^^^ to be aelf-averaging aa JV — k so. To date, a proof of thia property haa not been 
given, even though the adf-averaging of the analogoua quantity in the more general matching 
problem (where linka P — P and Q — Q arc allowed aa well) can be ahown at all d in caaentially 
the aame way aa for the TSP, fiollowing argumeuta developed by Steele [V]. Ear d = 1, it ia 
in fact known that aelf-averaging faili in the BM. Eor large d, however, let ua aaaume that 
Zbii/W'^"^^' doea converge to some /JHii('i) i^ the large N limit. 
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Wc aliAll Diow dmvc a boiLiul for the EncHdcan TSP coniitaiit /Se(<^ ^ tcmia of (iEui'^)' 
Caasadici K difijoiut sicta Si,. ^S;^ , togcthicr fbnmng a. Ijirgc sict S = Si U • " U S;^, ,iiid let 
csLch. Sj coatflia N r;mdiam pdnta in the tf-diiuciifiioaal unit hypcumhc. Confitnict the K 
iiniiiiTiniiTin matdiinga Si — Sa, Sa — Sa, - , Sk--! — S^r and Sk- — Si. Staitin^ At any point in 
Si, gencxAtc n loop (a clofiod path) in S by fcllowin^ the matdiin^ Si — Ss, Ss — S3, . until 
the path rc tninA to i\& fitartiitg point. The fiet of all fmch difitinct loop^ fli , . . . ,fljhr (M < N) 
\a then eqnii^knt to the siet S, and fiutthermiOite the anm of the loop len^tha ia oqnal to the 
finni of all nummnm matching eoatfi (Note that defined as 

Now, couiider the optimnm TSP tonr thiioi^halL the pointa of Si. Conatmet a giant doaed 
path viiiitiivg every podnt in S at least once, by snbstitntiqg into this TSP tonr the loops 
rii,...,njhr in place of their starting pointi in Si. Using standard techniques [G], wc can 
constmct from this path of length (iif -|- 1 )JV a shorter closed path of length KN which visits 
every point in S exactly once. For the Endidean TSP tonr lei^th £e, we then obtain the 
inequality 

K 

1=1 

If S consists of random points chosen independently and nnifomaly in the unit hypercnbe, 
then averaging over all configurations, dividiDg by N^~^^^ and taJdng the Unut N — i' ao, we 
find 

K ^ ^^.iSE (d) <(iE(d) + KfiEia (di . (D 2) 

Letting K = d^ this gives in the large d limit iii:{d) < /3Bu(d). Based on a nalogies with 
other combinatorial optimization ptoblems [l^], /Shu('^ expected to scale as •\/d/2mi when 
d — I- ao. In that case, (i^id) too must satisfy the E(ertsimas-van Ryzin conjecture (IG). 
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— Chapter overview — 

Thus far, we have used the Euchdean stochastic tsp as the basis for discussion. We have 
studied the behavior of the optimal tour length Le{N, d) at large N using numerical techniques. 
We have seen how, when N —>■ oo, LE{N,d)/N^'~^/'^ converges to its instance-independent limit 
f}E{d). We have then taken advantage of the random link approximation to obtain PRL{d), a 
theoretical prediction for f3E{d). This approximation involves substituting a different model for 
the Euclidean model: the random link TSP, in which the independent random variables are not 
the positions of cities but rather the lengths separating pairs of cities. 

The random link TSP was developed by theoreticians in search of an analytically tractable 
version of the problem. A major breakthrough occurred with the use of the cavity method, by 
Mezard & Parisi (19866) and later Krauth & Mezard (1989), to provide an analytical 
prediction for the (random link) quantity PuLid). In this chapter, we discuss the analytical 
approach afforded by the cavity method. The approach relies on certain assumptions which, 
while plausible, have thus far been tested numerically only for a uniform distribution of link 
lengths — corresponding to the d = 1 case. We therefore perform numerical simulations to 
test the cavity predictions more broadly, and find convincing evidence that these predictions are 
correct for all d. 

The strategy for attacking the random link TSP analytically (Orland, 1985) has been to 
map the problem onto a model with m-component spins; taking the limit m — > leads to 
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a representation of a self-avoiding walk. This method, originally developed in the context of 
polymer theory (De Gennes, 1972), then allows us to recover the partition function for the 
TSP in terms of a model that is more easily solvable. The first success in solving it was due to 
Mezard & Parisi (1986a), by means of the "replica" method developed in spin glasses (see 
Mezard, Parisi & Virasoro, 1987). The idea behind the replica method is a simple one. 
In order to compute the free energy F = —TlnZ averaged over the ensemble, we imagine n 
independent copies (replicas) of the system, and then write 

(InZ) = lim ^^^^ -. (III.l) 

Making the assumption that all n replicas are interchangeable, i.e., they all reach the same 
equilibrium state (this ergodicity hypothesis is known as replica symmetry), Mezard and Parisi 
obtained a high-temperature expansion for PRL{d =1). 

Unfortunately, the replica method allowed estimating the solution at T = (the global opti- 
mum) only by a numerical extrapolation from higher temperatures. Mezard and Parisi therefore 
instead tried using the cavity approach, a mean-field method (also assuming replica symmetry) 
whereby one adds an additional spin to the system and calculates the magnetization at this site 
in the absence of correlations between the other spins. This mean-field assumption is generally 
held to be correct in the limit of a large number of spins. Krauth & Mezard (1989) then 
completed the resulting calculation, finding an integral equation for (3RL{d). Cere, Boutet 
DE MONVEL, BOHIGAS, MARTIN & Percus (1997) solved the integral equation numerically for 
small values of d, as well as in a large d expansion: 



2-ln2-27 



(m.2) 



where 7 is Euler's constant. 

The problematic assumption used here, however, is that of replica symmetry. In the spin 
glass model of Sherrington & Kirkpatrick (1975), where similar methods were first tried, 
it was found that the replica symmetry assumption led to a ground-state energy that was off 
by about 5% (Mezard, Parisi & Virasoro, 1987, pp. 13-14). It was therefore necessary 
to resort to a replica symmetry breaking (rsb) scheme in order to obtain more realistic results. 
This involved defining an "overlap" Qab between two states a and b, so as to provide a measure of 
"distance" between them in state space. For the TSP, for example, Kirkpatrick & Toulouse 
(1985) have defined Qab as the proportion of links in common between two tours a and b. RSB 
implies, among other things, that q is not a self-averaging quantity, i.e., in the large N limit the 
distribution P{q) docs not simply approach a delta function. 

In the case of the random link tsp, Krauth & Mezard (1989) performed a number of 
d = 1 numerical checks indicating that, unlike for the spin glass, the replica symmetric solution 
given by the cavity method is correct. They solved the cavity equation numerically to obtain 
Prl{^) = 1-0208, finding it in close agreement with the earlier direct simulations of Kirkpatrick 
and Toulouse — whose results suggested /3rl(1) ~ 1.045 — as well as with Mezard and Parisi's 
extrapolated vahio of — 1.04 ± 0.015 from the replica method. Recent siniTilations by 
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Johnson, McGeoch & Rothberg (1996) give /3rl{1) = 1.0209 ± 0.0002 (see Appendix E), 
providing excellent confirmation of the theoretical predictions. Krauth and Mezard also per- 
formed simulations of their own, and compared the cavity and simulated values for the proba- 
bility distribution Vdil) of link lengths / in the optimal tour, finding close agreement. 

More direct evidence of the lack of RSB comes from the analysis by SOURLAS (1986), who 
considered the low-tcmpcrature statistical mechanics of the d = 1 random link TSP. Performing 
numerical simulations at temperatures as low as T = 0.85, he found that the distribution for the 
overlap function, as defined above, does indeed approach a delta function at large N. This sort 
of behavior is inconsistent with RSB, and suggests that the replica symmetric solution used in 
the cavity method is indeed exact. As Sourlas noted, however, obtaining good enough statistics 
for simulations is extremely difficult for lower T, where fluctuations in q become much larger. 

We wish to extend these numerical confirmations of the cavity method to d > 1, in order 
to determine its validity for choices of d relevant to the Euclidean tsp. Since we are interested 
in T = properties, direct analysis of the overlap function does not appear within reach. 
However, the sort of tests proposed by Krauth and Mezard is easily applicable to cases of 
interest such as d = 2, and another method can be used to test the 0{l/d) coefficient in the 
large d expansion (III. 2) for l3RL{d)- (The leading order term has already been shown exact 
by Vannimenus & Mezard (1984).) This latter method involves defining a "renormalized" 
random link TSP without the parameter d, having tour lengths whose leading order behavior at 
large N goes as L{N) ~ jiN. The quantity may be expressed in terms of the 0{l/d) coefficient 
in (III. 2); by performing numerical simulations on this renormalized model, we manage to verify 
the subleading behavior of the cavity method's /3_rl(c?) prediction. 

Our simulation methods are summarized in Appendix D and follow those of Chapter II: we 
use heuristic algorithms at values of A'^ up to A'^ = 100, where both statistical and systematic 
errors can be well controlled. We extrapolate to the large N limit by fitting to the expected 
finite size scaling law. The resulting fits are excellent, and confirm the cavity predictions both 
of Prl{^) and of the subleading term in the large d (3RL{d) expansion to within well under 1%. 
(This may be compared with the 5% error of the replica symmetric solution in the Sherrington- 
Kirkpatrick spin glass model.) Taken together, then, this provides good arguments for the claim 
that the cavity method is exact for all d, and further reason for conjecturing that there is no 
RSB down to T = 0. 
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Abstract 

We study the random hnk TSP, where lengths kj between city i and city j {i < j) 
are taken to be independent, identically distributed random variables. We discuss 
a theoretical approach, the cavity method, that has been proposed for finding the 
optimal tour length over this random ensemble, given the assumption of replica 
symmetry. Using finite size scaling and a renormalized model, wc test the cavity 
predictions against the results of simulations, and find excellent agreement over a 
range of distributions. In doing so, we provide numerical evidence that the replica 
symmetric solution to this problem is the correct one. 
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1. Introduction 



— 1 — Introduction 

Over the past 15 years, the study of the travehng salesman problem (tsp) from the point of 
view of statistical physics has been gaining added currency, as theoreticians have improved 
their understanding of the links between combinatorial optimization and disordered systems. 
One particular breakthrough occurred with the idea, first formulated by Mezard & Parisi 
(19866) and later developed by Krauth & Mezard (1989), that a version of the stochastic 
TSP could be solved using an analytical method inspired from spin glasses. This method, known 
as the cavity method, is based on certain assumptions pertaining to the physical properties of 
the system, notably that of replica symmetry. Although in the case of the spin glass replica 
symmetry is known to be broken (Mezard, Parisi &: ViRASORO, 1987), for the tsp there 
are various grounds for suspecting that the assumption is valid (SOURLAS, 1986; Krauth & 
Mezard, 1989). 

The traveling salesman problem is as follows: given N sites (or "cities"), find the length 
of the shortest closed path ("tour") passing through all cities exactly once. In the stochastic 
TSP, the matrix of distances separating pairs of cities is drawn randomly from an ensemble. 
The ensemble that has received the most attention among theoreticians is the random link case, 
where the lengths lij between city i and city j {i < j) are taken to be independent random 
variables, all identically distributed according to some p{l). The idea of looking at this random 
link ensemble, rather than the more traditional "random point" ensemble with cities distributed 
uniformly in Euclidean space, originated with an attempt by Kirkpatrick to find some sort of 
TSP equivalent to the earlier Sherrington & Kirkpatrick (1975) model for spin glasses. 

The great advantage of the random link TSP over the (random point) Euclidean TSP is that 
with the former it is actually possible to make analytical progress on the problem. The cavity 
solution of Mezard & Parisi (19866) and Krauth & Mezard (1989) leads to a system of 
integral equations, which can be solved — numerically at least — to give the optimal tour length 
in the limit N ^ 00. 

In a previous paper (Cerf, Boutet de Monvel, Bohigas, Martin & Pergus, 1997), we 
chose the random link distribution p{l) to match that of the distribution of individual city-to-city 
distances in the Euclidean case, and used the random link TSP as a random link approximation 
to the Euclidean tsp. The approximation might seem somewhat crude, since it neglects all 
correlations between Euclidean distances, such as the triangle inequality. Nevertheless, it gives 
remarkably good results. For the d-dimensional Euclidean case in d = 2 and d = 3, a numerical 
solution of the (random link) cavity equations predicts asymptotic N ^ 00 optimal lengths 
within 2% of the values obtained from direct (Euclidean) simulation. In the limit d ^ 00, this 
gap shows all signs of disappearing. The random link problem is thus more closely related to 
the Euclidean problem than most researchers have suspected. 

The random link TSP is also interesting in itself, however, particularly because little numerical 
work has accompanied the analytical progress made. This is all the more important given the 
questions surrounding the hypotheses adopted in the cavity method. In this paper we attempt 
to redress this imbalance, providing first of all a numerical study of the finite size scaling of the 
random link optimal tour length, and second of all. empirical arguments suggesting that the 



55 



Chap. Ill: The random link TSP and its analytical solution 



cavity solution is in fact the correct one. 



— 2 — Background and the cavity method 

In an attempt to apply tools from statistical mechanics to optimization problems, KiRKPATRiCK 
& Toulouse (1985) introduced a particularly simple case of the random link TSP. The distri- 
bution of link lengths lij was taken to be uniform, so that p{l) is constant over a fixed interval. 
In light of the random link approximation, one may think of this as corresponding, at large 
A^, to the 1-D Euclidean case. (When cities are randomly and uniformly distributed on a line 
segment, the distribution of lengths between pairs is uniform.) Although the 1-D Euclidean case 
is trivial — particularly if we adopt periodic boundary conditions, in which case the optimal 
tour length is simply the length of the line segment — the corresponding random link problem 
is far from trivial. 

The simulations performed by Kirkpatrick and Toulouse suggested a random link optimal 
tour length value of Lrl ~ 1.045 in the N ^ 00 hmit.^ Mezard & Parisi (1986a) attempted 
to improve both upon this estimate and upon the theory by investigating the random link TSP 
using replica techniques often employed in spin glass problems. This approach allowed them to 
obtain, via a saddle point approximation, many orders of the high-temperature expansion for 
the internal energy. They then extrapolated down to zero temperature — corresponding to the 
global TSP optimum — finding Lrl = 1-04 it 0.015. This analysis, like that of Kirkpatrick and 
Toulouse, was carried out only for the case of constant p{l). 

Given the difficulties of pushing the replica method further, Mezard and Parisi then tried a 
different but related approach known as the cavity method (Mezard & Parisi, 19866). This 
uses a mean-field approximation which, in the case of spin glasses, gives the same result as the 
replica method in the thermodynamic limit (iV — > 00). Let us sketch what is involved in the 
cavity method, not least so that we may enumerate clearly the assumptions made. 

Both the replica and the cavity approaches involve mapping the TSP onto an m-component 
spin system, writing down the partition function at temperature T, and then taking the limit 
m — 0. More explicity, consider N spins Si, i = 1, . . . ,N (corresponding to the TV cities), where 
each spin Sj has m components S^, a = 1, . . . ,m, and where Sj • Sj = m for all i. The partition 
function is defined, in terms of a parameter lo, as 

Z = j dn{S} exp(u; ^ Rij Sj • S^) (III.3) 



i<3 



j dii{S] 



1 + u; ^ i?i,(Si • S,) + ^ ^ Rij RkiiSi ■ Sj){Sk ■ Si) + 



2! 

k<l 



(III.4) 



^Here we work in units where the hne segment is taken to have unit length, and in order to match the 
normahzed 1-D Euchdean distribution, we let p{l) = 2 on (0,1/2). The 1-D Euclidean value, for comparison, 
would thus be Le — 1. Kirkpatrick and Toulouse, among others, choose instead p{l) = 1 on (0, 1), contributing 

an additional factor of 2 in Lbl wliic'li wo omit when quoting' their rosnlts. 
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where J dn{S} denotes an integral over all possible spin values (this simply corresponds to the 
surface of an m-dimensional sphere), and Rij is related to the length lij between city i and 
city j as Rij = e~^^^''^'-^^'^ . Now we employ a classic diagrammatic argument: let each spin 
product (Sp • Sq) appearing in the series be represented by an edge in a graph. The first-order 
terms (uj) will consist of one-edge diagrams, the second-order terms {uj'^) will consist of two- 
edge diagrams, and so on (see Figure III-l). What then happens when we integrate over all 
spin configurations? If there is a spin Sp which occurs only once in a given diagram, i.e., it 
is an endpoint, the spherical symmetry of Sp will cause the whole expression to vanish. We 
will therefore be left only with "closed" diagrams, where there is at least one loop. It may 
furthermore be shown that in performing the integration, any one of these closed diagrams will 
contribute a factor m for every loop present in the diagram (De Gennes, 1972; Orland, 1985). 
If we then consider {Z — l)/rn and take the limit m 0, it is clear that only diagrams with a 
single loop will remain. Furthermore, since any closed diagram with more than N links must 
necessarily contain more than one loop, only diagrams up to order will remain. Finally, take 
the limit u — ^ 00. The term that will then dominate in (111.4) is the order 00^ term which, being 
a single loop diagram, represents precisely a closed tour passing through all sites. We may 
write it without the combinatorial factor A'^! by expressing it as a sum over ordered pairs in the 
tour, and we thus find: 



Z — \ ^ 

W ~ / V ^ii2 -^'his ' ' ' ^iN-iiN ^iiiN (III. 5) 

m—^U TnU! 

(jj-^oo W-link single loops 

(il, 22, •••!««) 

A''-city tours 

where L is the total tour length. We thus obtain exactly the partition function for the traveling 
salesman problem, with the correct canonical ensemble Boltzmann weights, using the tour length 
as the energy to be minimized (up to a factor N^/'^, necessary for the energy to be extensive). 



Figure III— 1: A graphical representation of a 1-edge open diagram; a 2-edge open diagram; an 8-edge 
closed diagram with two loops; and a 12-edge closed diagram with one loop. 
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The cavity method now consists of adding an additional (A'^ + l)th spin to the system in the 
presence of an infinitesimally small external magnetic field oriented along component 1, leading 
to a spontaneous magnetization (Sl). The magnetization {Slf_^i) is then expressed in terms of 
what the other {SD^s would be in the absence of site A'^ + 1, hence the notion of a "cavity". 
In order to do this, two important assumptions are made. First of all, wc assume that there is 
one single equilibrium state. That is to say, at all temperatures a unique thermodynamic limit 
exists. A similar assumption of ergodicity is made to obtain the replica solution given earlier. 
The assumption is known as replica symmetry: if a large number of "copies" of the system are 
allowed to develop, they all reach the same equilibrium state. The notion of replica symmetry 
is central to the study of disordered systems. Replica symmetry is known to be broken in spin 
glasses, though there is no reason why this should necessarily be so in other NP-hard problems. 
Showing that the cavity solution correctly predicts macroscopic quantities for the random link 
TSP suggests that indeed the tsp exhibits replica symmetry. 

The second important assumption on which the cavity method is based involves neglecting 
the correlations between spins Si,i = 1,...,N, when calculating {Slf_^i). The assumption is 
that at large A'", the contributions due to these correlations are 0{1/N) or less, so that when 
N ^ oo this sort of mean- field approximation is justified. Proceeding in this way, one obtains 
from the partition function a recursion relation for {Sl^_^^) in terms of the other {Sj)^s: 



^l<i<j<N RN+l,i{Sl)RN+l,j{Sj) 



Furthermore, it may be seen from (III. 4) that the diagrams contributing to the thermal average 
(Sp • Sq) are those that have link p — q occupied; (Sp • S^) is then the occupation number Upq of 
the link. Since the spontaneous nia,gii6tizcitioii is cilong component 1, Tipq is simply (SpSg). From 
the partition function, a recursion relation for nN+i,i in terms of the (5'/)'s may be obtained 
much as in (III. 7): 

nM-,i,i = RN^iASU,){Sl)^^^fi^^^^. (III.8) 

Now consider the effect of our second (mean-field) assumption over the ensemble of instances 
(distribution over the disorder). As far as (HI. 7) is concerned, we may treat the magnetizations 
(Sj) as independent identically distributed random variables. Requiring {Slf_^_-^) to have the 
same distribution as the others imposes, for a given link length distribution p{l), a unique self- 
consistent probability distribution of the magnetizations. Via (III. 8), then, the distribution 
Vd{l) of link lengths I in the optimal tour (at A^ oo) may be found. Krauth & Mezard 
(1989) carried out this calculation in the T — limit, for p{l) corresponding to that of the 
d-dimensional Euclidean case, namely 
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They found 



where Hd{x) is the solution to the integral equation 

Prom Vd{l), one may obtain the mean Hnk length in the tour, and thus the total length of the 
tour: 

r+00 

LRL{N,d) = N lVd{l)dl (III.13) 
^0 

d f+°° 

= TV^-Vrf^ / Hdix)[l + Hdix)]e-^''^''Ux (III.14) 
2 J -00 

for large N. At d = 1, Krauth and Mezard solved these equations numerically, obtaining Lrl = 
1.0208. It is difficult to compare this with Kirkpatrick's value of 1.045 from direct simulations (as 
no error estimate exists for the latter quantity), however Johnson, McGeoch & ROTHBERG 
(1996) have recently obtained the numerical result Lrl{N = 10,000) = 1.0211 it 0.0003, giving 
strong credence to the cavity value. Krauth and Mezard also performed a numerical study of 
'Pi{l). They found the cavity predictions to be in good agreement with what they found in 
their own direct simulations. Further numerical evidence supporting the assumption of replica 
symmetry was found in an analysis, by SOURLAS (1986), of the low temperature statistical 
mechanics of the system. Thus, for the lij distribution at d = 1, there is good reason to believe 
that the cavity assumptions are valid and that the resulting predictions are exact at large N . 

The cavity solution was extended to higher dimensions by the present authors (Fergus & 
Martin, 1996; Cerf, Boutet de Monvel, Bohigas, Martin & Fergus, 1997), and a large 
d power series solution was given for the asymptotic value PRiid) = limjv^oo Lrl^N^ d)/N^~^^'^: 



2 - In 2 - 27 _ / 1 



(III.15) 



where 7 represents Euler's constant (7 ~ 0.57722). It would be nice, then, to have evidence 
that the cavity method is exact for all d and not just at = 1. While certain arguments 
have been advanced in favor of this claim (Mezard, Farisi & ViRASORO, 1987; SouRLAS, 
1986), they have yet to be backed up by numerical evidence — as they have been, for instance, 
in a related combinatorial optimization problem known as the matching problem (Brunetti, 
Krauth, Mezard & Farisi, 1991; Boutet de Monvel, 1996). We now turn to this task, 
considering first the d = 2 case, and then a "renormalized" random link model enabling us to 
verify numerically the 0(l/d) coefficient in (III. 15). 
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— 3 — Numerical analysis: d = 2 case 

We have implicitly been making the assumption so far, via our notation, that as A'^ ^ oo 
the random variable LRL{N,d)/N^^^/'^ approaches a unique value fSmid) with probability 1. 
This is a property known as self-averaging. The analogous property has been shown for the 
Euclidean TSP at all dimensions (Beardwood, Halton & Hammersley, 1959). For the 
random link TSP, however, the only case where a proof of self-averaging is known is in the c/ — ^ 00 
limit, where a converging upper and lower bound give in fact the exact result (Vannimenus & 
Mezard, 1984): 



i + o 



(III. 16) 



Note, incidentally, that this already shows that the cavity solution (III. 15) is correct in the 
infinite dimensional limit. 

For finite d, however, it has not been shown analytically that Pmid) even exists. To some 
extent, the difficulty in proving this can be traced to the non-satisfaction of the triangle inequal- 
ity. The reader acquainted with the proof by Beardwood, et al, may see that the ideas used 
there are not applicable to the random link case; for instance, combining good subtours using 



0.2 




Rescaled tour length 

Figure III-2: Distribution of 2-D random link rescaled tour length (L/fx — {Lrl))/"/N for increasing 
values of N. Plus signs show N = 12 (100,000 instances used), squares show = 17 (100,000 instances 
used), diamonds show A'' = 30 (4,000 instances used), and dots show N = 100 (1,200 instances used). 
Solid lines represent Gaussian fits for each value of N plotted. 
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simple insertions will not lead to near-optimal global tours, making the problem particularly 
challenging. Let us therefore examine the distribution of d = 2 optimum tour lengths using 
numerical simulations, in order to give empirical support for the assertion that the AT — > oo 
limit is well-defined. 

In Figure III-2, we see that the Lrl{N,2)/\/N distribution indeed becomes increasingly 
sharply peaked for increasing A'^. Furthermore the variance of Lf{L{N,2) remains relatively 
constant in N (see Table III-l), indicating that the width a for the distribution shown in the 
figure decreases as l/\/N, strongly suggesting a Gaussian distribution. Similar results were 
found in an earlier study of ours concerning the Euclidean TSP (Cerf, Boutet de Monvel, 
BoHiGAS, Martin & Fergus, 1997) (albeit in that case with a approximately half of its 
random link value). This is precisely the sort of behavior one would expect were the central 
limit theorem to be applicable. 

Table III-l: Variance of the non-rescaled optimum tour length Ljii^{N, 2) with increasing A'^. 



N 


(7^ 


# instances used 


12 


0.3200 


100,000 


17 


0.3578 


100,000 


30 


0.3492 


4,000 


100 


0.3490 


1,200 



It is worth noting that self-averaging in a different quantity, the Farisi overlap q (Mezard, 
Farisi, Sourlas, Toulouse & Virasoro, 1984; Mezard, Farisi & Virasoro, 1987), would 
in fact be a direct indication of replica symmetry. The overlap Qah between two states (tours) 
a and b may be defined in the TSP, following Kirkpatrick and Toulouse, as the fraction of links 
that are common to tours a and h. If we then consider suboptimal tours produced by a finite- 
temperature algorithm, and measure the overlaps between these tours, replica symmetry requires 
that the overlap distribution become more and more sharply peaked as A'^ ^ 00. Results by 
Sourlas (1986) confirm that this is so for d = 1, at temperatures going down to T = 0.85. 
Unfortunately, such simulations at lower temperatures do not appear feasible: Sourlas' results 
suggest that fluctuations in q become uncontrollable as T — 0. As we are interested primarily 
in the T = case, we do not consider q in the present analysis, resorting instead to direct tests 
of the cavity predictions. 

The algorithmic procedures we use for simulations are identical to those we used in our 
Euclidean study; for details, the interested reader is referred to that article. Let us, however, 
note here that our optimization procedure involves using the lk and CLO local search heuristic 
algorithms (LiN &; Kernighan, 1973; Martin & Otto, 1996) where for each instance of the 
ensemble we run the heuristic over multiple random starts. LK is used for smaller values of N 
{N < 17) and GLO, a more sophisticated method combining LK optimization with random jumps, 
for larger values of N {N = 30 and N = 100). Fractitioners may note at this point that given 
the problems we have mentioned in obtaining near-optimal tours from near-optimal subtours in 
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0.726 




0.02 0.04 0.06 0.08 0.1 
l/N 

Figure III-3: Finite size scaling of d = 2 optimum. Best fit (x^ ~ 4.46) is given by: {Lhl{N, 2)) /N^/^ = 
0.7243(1 + 0.0322/A^ - 1.886/N^). Error bars sliow one standard deviation (statistical error). 



the random link TSP, local search methods would be expected to perform very poorly leading to 
relative excess lengths which diverge as ^ oo and thus next to useless for a large N analysis. 
In the case of the lk heuristic, interestingly enough, the divergence appears to be no worse than 
logarithmic in (JOHNSON &: McGeoch, 1997), presumably because of the "variable depth" 
search that lk conducts. There is, nevertheless, a certain probability that even over the course 
of multiple random starts, our heuristics will not find the true optimum of an instance. We 
estimate the associated systematic bias using a number of test instances, and adjust the number 
of random starts to keep this bias at least an order of magnitude below other sources of error 
discussed below. (At its maximum — occurring in the = 100 case — the systematic bias is 
estimated as under 1 part in 20,000.) 

Following this numerical method, let us now consider the large limit of LfiL{N, 2). In the 
Euclidean case, it has been observed that the finite size scaling law can be written in terms of 
a power series in 1 /N. We expect this same behavior in the random link case, namely that the 
ensemble average {Ljil{N ,2)) satisfies 

{LRL{N,d))=f3RL{d)N^-^l'' 

In order to obtain {Lrl{N, 2)) at a given value of A^ from simulations, we average over a large 



1 + 



m 

N 



+ 



(III.17) 
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number of instances to reduce the statistical error arising from instance-to-instance fluctuations. 
Figure III-3 shows the results of this, with accompanying error bars, fitted to the expected finite 
size scaling law. The fit is a good one (x^ = 4.46 for 5 degrees of freedom) and gives an extrap- 
olated value of pBLi^) = 0.7243 it 0.0004. This is in superb agreement with the cavity result 
of 0.7251. The relative discrepancy between the cavity prediction and our numerical estimate 
is approximately 0.1%, which in consistent with our statistical error bars; by comparison, the 
error in the replica symmetric solution to the Sherrington-Kirkpatrick spin glass ground state 
energy is estimated to be of the order of 5% (Mezard, Parisi & ViRASORO, 1987, pp. 13-14). 

Another quantity that was considered in the d = 1 numerical study of Krauth & Mezard 
(1989) is the optimal tour link length distribution Pail) given in (III. 10). Let us consider 
^2(0) ^^"^ following their example, let us look specifically at the integrated distribution Id{l) = 
J^Vdil) dl. The cavity result for Idil) can, like fimid)-, be computed numerically to arbitrary 
precision. In Figure III-4 we compare this with the results of direct simulations, at increasing 
values of N . The improving agreement for increasing N (already within 2% at TV" = 100) strongly 
suggests that the cavity solution gives the exact N ^ 00 result. 

Finally, it may be of interest to consider one further quantity in the d = 2 random link 
simulations: the frequencies of "neighborhoods" used in the optimal tour, that is, the propor- 
tion of links connecting nearest neighbors, 2nd-nearest neighbors, etc. While there is no cavity 
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Figure III-4: Integrated probability distribution of link lengths in the optimal tour, for d = 2, using 
rescaled length V = iVN. Plus signs represent N = 12 simulation results, dots represent N = 100 
simulation results, and solid line represents cavity prediction. 
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Figure III-5: Frequencies with which fcth- nearest neighbors are used in optimal 2-D random hnk tours. 

Plus signs show values for = 12, squares for N = 17, diamonds for N = 30, and dots for A'' = 100. 
Best exponential fit (appearing linear on log plot) is shown for N = 100 data. 

prediction for this quantity, SOURLAS (1986) has noted that in practice in the d = 1 case, the 
frequency falls off rapidly with increasing neighborhood — suggesting that optimization heuris- 
tics could be improved by preferentially choosing links between very near neighbors. Results of 
our simulations, shown in Figure III-5, suggest that at d = 2 the decrease is in fact very close 
to exponential. We have no theoretical explanation for this behavior at present, and consider it 
a fascinating open question. 

— 4 — Numerical analysis: renormalized model 

Let us now consider a different sort of random link TSP, proposed in Cerf, Boutet de Mon- 
VEL, BOHIGAS, Martin & Fergus (1997), allowing us to test numerically the l/d coefficient 
in the cavity result (III.15). Define {Di{N,d)) to be the distance between a city and its nearest 
neighbor, averaged over all cities in the instance and over all instances in the ensemble.^ For 
large d, it may be shown that 



lim N^/'^{Di{N,d)) 




) 



(III. 18) 



^N()t<> that (Di_(N. d)) itsc;lf docs not involve the notion of optimal tonrs. or tours of any sort for that matter. 
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where 7 is Euler's constant. It is not surprising that this quantity is reminiscent of (III. 16), 
since N^^'''{Di{N,d)) represents precisely a lower bound on PRiid). 

In the renormalized random link model, we define a new link length Xij = d[lij — 
{Di(AL,d))]/{Di{N,d)), where the kj have the usual distribution corresponding to d dimen- 
sions. The Xij are "lengths" only in the loosest sense, as they can be both positive and negative. 
The optimal tour in the Xij model will, however, follow the same "path" as the optimal tour in 
the associated lij model since the transformation is linear. Its length Lx{N,d) will simply be 
given in terms of Li{N, d) by: 



{Di{N,d)) 



Li{N,d) = N{DiiN,d)) 



1 + 



LAN,d) 
dN 



By definition, then, 

pRLid) 



lim iV^/'^(L>i(iV,d)), or at large d, 

N^oo 



1 



7 



4 + 



d 



1 



d2 



lim 



1 + 



dN 



(III.19) 
(III.20) 

(III.21) 
(111.22) 



If PRL{d) is to be well-defined, then there must exist a value fj,{d) such that 
lim^v^oo Lx{N, d)/N = n{d). 

What will be the distribution of lengths p{x) corresponding to p(/)? Prom the definition of 
p{l) and that of x, 



p{x) 



, and substituting for /, 



r{d/2 + 1) d 

tt''/- / r\d-l 



In the limit AT — > 00, this gives 

, ^. . j / ^ yV'i 



p{x) 



N- 



-'('-F( 



d-l 



1 + 



by Stirling's formula 



1 + 



(III.23) 
(III.24) 

(III.25) 
(III.26) 
(III.27) 



At large d, we sec that p{x) is to leading order independent of d; the same must then be true for 
Lx{N,d), and hence for p{d), which we now write as po[l + 0{l/d)]. From (III. 22), we obtain 



1 + 



^0-7 

d 



(III.28) 



Finally, we may perfectly well define p{x) in the d — > 00 limit, in which case (III. 27) gives 

p(x) = A^^"^ oxp(,T — 7). Doing so will give us a renormalized model with lini/v^oc Lj-{N)/N = fiQ. 
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By performing direct simulations on a random link model with this exponential distribution, 
we may find the value of //q numerically, and thus from (III. 28) the 1/d coefEcient for the 
d-dimensional (non-renormalized) case. 

Figures III-6 and III-7 show the numerical results for the renormalized model. In Figure 

III-6, we see that just as in the d = 2 case, the distribution of the optimal tour length becomes 
sharply peaked at large N and the asymptotic limit jiQ is well-defined. Via (III. 22), this provides 
very good reason for believing that (3RL{d) is well-defined for all d, and self-averaging holds for 
the random link TSP in general. In Figure III-7, we show the finite size scaling of {Lx{N)/N) . 
The fit is again quite satisfactory (with x^=5.23 for 5 degrees of freedom), giving the asymptotic 
result /zq = 0.7300 ± 0.0010. The simulated value for the 1/d coefficient in Pniid) is then 
HQ — J = 0.1528 lb 0.0010, in excellent agreement (error under 0.3%) with the cavity value 
2 - In 2 - 27 0.1524 given in (III. 15). 

Also as in the d = 2 case, let us briefly consider the frequencies of fcth-nearest neighbors 
used in optimal tours. This is given in Figure III-8. Even though the exponential fit is not 
as excellent as in the d = 2 case, it is still striking here. We have seen that the renormalized 
random link model is analogous to the d ^ 00 limit of the standard random model, as far as 
tour paths are concerned. These fcth-neighbor frequency results then give empirical evidence 



0.2 




Rescaled tour length 

Figure III-6: Distribution of renormalized random link rescaled tour length [L^ — {L^) ) /N for increasing 
values of N. Plus signs show = 12 (100,000 instances used), squares show iV = 17 (100,000 instances 
used), diamonds show A'' = 30 (4,000 instances used), and dots show N = 100 (1,200 instances used). 
Solid lines represent Gaussian fits for each value of N plotted. 
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0.742 
0.738 
0.734 

1 

0.73 
0.726 
0.722 

0.02 0.04 0.06 0.08 0.1 

Figure III-7: Finite size scaling of "renormalized model" optimum. Best fit (x^ = 5.23) is given by: 
{L^{N))/N = 0.7300(1 + 0.3575/iV - 2.791/iV^). Error bars show one standard deviation (statistical 
error) . 

that even at infinite dimension, the "typical" k used remains bounded. 

— 5 — Conclusion 

The random link TSP has interested theoreticians primarily because of its analytical tractability, 
allowing presumably exact results that are not possible in the more traditional Euclidean tsp. 
Other than in the d = \ case, however, it has attracted little attention from the angle of 
numerics. In this paper we have provided a numerical study of the random link TSP that 
was lacking up to this point, and which addresses some long-unanswered questions. Through 
simulations, we have tested the validity of the theoretical cavity method predictions. While in 
other disordered systems, such as spin glasses, the replica symmetric solution gives values of 
macroscopic quantities that are inexact (typically by at least 5%), in the random link tsp it 
shows all signs of being exact. We have studied various macroscopic quantities at d = 2 and 
found that the numerical results confirm the cavity predictions to within 0.1%. Furthermore, we 
have confirmed, by way of a renormalized random link model, that the analytical cavity solution 
gives a large d expansion for the optimal tour length whose 1 /d coefficient is correct to within well 




67 



Chap. Ill: The random link TSP and its analytical solution 




■I^Q"-" I .... I .... I .... I .... I .... I .... I .... I .... I .... I 

123456789 10 
Neighborhood (k) 

Figure 111—8: Frequencies with which fcth- nearest neighbors are used in optimal renormaUzed random 
link tours. Plus signs show values for A'^ = 12, squares for A'' = 17, diamonds for N = 30, and dots for 
A'' = 100. Best exponential fit (appearing linear on log plot) is shown for N = 100 data. 

under 1%. The excellent agreement found at cZ = 1 (Krauth & Mezard, 1989; SOURLAS, 1986), 

d = 2, and to 0{l/d) at large d, then suggest strongly that the cavity predictions are exact. 
This, finally, provides indirect evidence that the assumption of replica symmetry — on which 
the cavity approach is based — is indeed justified. 
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— Chapter overview — 



In this final chapter wc turn to a problem that does not involve the TSP as such, but is 
closely inspired by our tsp work and mirrors much of the language we have employed so far. 
We have seen in Chapter II how the N ^ oo Euclidean optimal tour length, where N repre- 
sents the number of cities, could be extrapolated from reasonably small sized instances via an 
understanding of the finite size scaling behavior. We were thus able to obtain precise numerical 
estimates for this asymptotic value without undue computational effort. 

In order to understand the finite size scaling, we decided to set aside temporarily the notion 
of optimal tours — or, for that matter, tours of any sort. Instead we simply looked at cities 
placed randomly and independently, with a uniform distribution, and considered the distances 
between nearest neighboring cities, second-nearest neighboring cities, and so on. What we found 
was an unexpected sort of universality. If one considers the ensemble average {Dk{N)) of the 
distances between /cth-nearest neighbors, the large N scaling law of this quantity is independent 
of k. Putting it another way, the A^-dependence and the fc-dependence of {Dk(N)) separate; in 
d dimensions we are left with, up to corrections exponentially small in N, 

^^'^^^^ iW~ r(iv + i/d) ■ ^^^-^^ 



Chap. IV: Universalities in nearest neighbor distances 



We now wish to investigate this universahty further. How, first of ah, does the property 
depend on the nature of the physical space used? The universality only applies in flat (i.e., 
Euclidean) spaces, but leads us to notice another interesting and more subtle property on curved 
surfaces. Let us look at the large behavior of {Dk(N)) (averaged over the whole surface, if 
the surface is not homogeneous). We find that {Dk(N)) scales as a power of times a series 
in whose 0{1/N) term, while now explicitly depending on k, has the form of a topological 
invariant: (x(2fc + 1) — 9)/24, where x is the Euler characteristic. To 0(1/N), then, it is the 
topology of the surface and not its detailed shape that plays a role in the finite size scaling 
behavior. 

The series describing {Dk{N)) may, for a general manifold, be obtained directly from 
the area A(l) of a geodesic disc of radius / on the manifold. A geodesic disc about some point 
X is simply defined as the locus of points whose distance from x along the manifold is smaller 
than the radius. For a flat surface, of course, A{1) = i^f'. For a curved surface, A(l) may be 
written as a series expansion in l\ higher-order terms in this expansion will then correspond to 
higher-order terms in the 1/N series. We must therefore find the series expansion of A{1), whose 
coefficients involve gradients of the surface's Gaussian curvature K at the point x about which 
the expansion is being performed. Unfortunately, we know of no general expression for these 
coefficients. Calculating them is a somewhat laborious procedure, which we have carried out up 
to the 0{l^) term in A{1): 

A{1) = T,f-llI-K + —(2K''-3V^K) 
^ ^ 12 720 ^ ' 

-— — — (8i^3 - 3[10(VK)2 + UKV^K - BV^K]) + 0(l^°). (IV.2) 
161280 

From this, terms through 0{l/N^) can be found in the 1/N series for {Dk{N)). We discuss 
these higher-order terms, although it appears that no further universalities exist beyond the 
topological invariant at 0{1/N). We then turn briefly to higher-dimensional manifolds — where 
A{1) is far less obvious — and examine the cases in which a topological invariant could exist 
there. Finally, we provide an interpretation of our results via a Regge calculus approach, where 
we consider sites distributed on a polyhedral surface rather than a smooth manifold; we find 
that there too, we recover the topological invariant. 

Let us alert the reader to some notational differences in this final chapter that we unfortu- 
nately cannot avoid. We find it clearer in the present context to speak of the distance from an 
arbitrary point x on the manifold to its /cth-nearest site, rather than of the distance between 
a site and its fcth- nearest neighbor. Df;{N) is thus a random variable defined everywhere in 
(continuous) space, and not only at the positions of one of the N sites. {Dk{N)) is its ensemble 
average, which for inhomogeneous surfaces will vary with x. Our previous notation amounted 
to restricting x to be the position of one of the cities, thus reducing the number of possible 
neighbors by 1. Our results for N cities, in our previous notation, are therefore equivalent to 
our results for A'^ — 1 cities, in our new notation. Note, also, that while the terms "cities" and 
"sites" are interchangeable, they are to be distinguished from the term "point", which is used 
here in its most precise geometric meaning to represent a,ny position x on the manifold. 
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Abstract 

Take N sites distributed randomly and uniformly on a closed surface. We express the 
expected distance {Dk{N)) from an arbitrary point on the surface to its /cth-nearest 
neighboring site, in terms of the surface area A{1) of a disc of radius / about that 
point. We then find two universalities. First, for a flat surface, where A{1) = ttZ^, 
the finite size scaling series giving corrections to the large N asymptotic behavior 
of {Di^{N)) is independent of k. Second, for a curved surface, the finite size scaling 
scries for the average J {Dk(N)) dfi over the surface is, to 0{1/N), a topological 
invariant. We discuss the case of higher dimensions (d > 2), and also interpret our 
results using Regge calculus. 



* Unite de recherche des Universites de Paris xi et Paris vi associee au CNRS 



1. Introduction 



— 1 — Introduction 

In many physics and computer science problems, a quantity of interest is the distance between 
neighboring sites in a space. One frequently wishes to calculate distances to a nearest neigh- 
bor, second-nearest neighbor, and more generally, A;th-nearest neighbor. Examples in computa- 
tional geometry and optimization abound, ranging from random packing of spheres to minimum 
spanning trees. Such problems also arise naturally in physics, both for interacting particle sys- 
tems such as liquids, and for cellular objects such as foams and random lattices (ITZYKSON & 
Drouffe, 1989). 

Here we consider the case of N sites placed randomly, with a uniform distribution, on a 
2-D surface of fixed area. Let the random variable Di.(N) represent the distance between a 
given point x and its fcth-nearest neighboring site. Its expectation value (Z)fc(A^)) taken over 
the ensemble of randomly placed sites — and in fact all moments {D^{N)) — then exhibit some 
surprising properties. 

When the surface is flat, {Dk{N)) at large N is described by its leading asymptotic behavior 
times a power series in where all orders of the series turn out to be independent of k. In 

the language of statistical physics, the finite size scaling series for {Dk{N)) is the same at all k 
to all orders in Geometrically, this universality is far from obvious. Furthermore, we find 

that the property is not restricted to two dimensions, and is equally valid for fiat spaces of any 
dimension. 

When the surface is curved, this ^-independence no longer holds, but one finds on the other 
hand another universality: the finite size scaling series, averaged over the entire surface, gives 
a topological invariant to 0{1/N). This term thus depends not on the detailed shape of the 
surface but only on the surface's genus. 

In this paper we explore these universalities. First, we express {Dk{N)) in terms of the 
area, A{1), of a disc of radius I on an arbitrary surface. We observe that in the special case 
where A(l) consists only of a power of I, the series for {Dk{N)) exhibits the universality 
in k, i.e., is independent of k. Second, we give the relation between A(l) and the Gaussian 
curvature of the surface, and find the leading correction in the power series for {Dk{N)). 
as functions of the curvature. This leads to the topological invariance at 0{1/N). We then 
discuss higher order terms, and the case of higher dimensions. Finally, we show that a Regge 
calculus approach provides a simple means of obtaining this topological invariance for the case 
of polyhedral (non-smooth) surfaces. 

— 2 — Preliminaries 

Take any point x, and consider P[Dk(N) = i]x, the probability density that the point's fcth- 
nearest neighboring site lies at a distance I from it. This is equal to the probability density of 
having k — 1 (out of N) sites within distance I, one site (out of N — k + 1) at distance I, and 
the remaining N — k sites beyond distance I. Let us choose units so that our surface has area 
1. Since sites are distributed uniformly over the surface, the probability of a single site lying 
within distance / is then simply the area ^(/)x of a disc of radiiis / about point x on the surface. 
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Dropping the argument x (in order to simpUfy the notation), we may then write 

giving the expectation value (first moment) 

{Dk{N)) = / P[DkiN) = l]ldl 
Jo 

^' r,r..nife-iri ,n^^N-k dA{l) 



{N-k)lik 



— l^^imr^ii-Aiir-'^^di. 



Under the variable transformation w = A(l), this may be written in terms of the inverse function 
A~^{w) as 

If A~^{w) admits the power series expansion in w: 

oo 

A~^{w) 2± w"' ^Cjw', for some 7 G [0, 1), (IV.3) 

j=0 



then 

ATI °° 

^'''^''^^ " (iV-fc)!(fc-i)! E X ^'^^'^^-'(1 - ^)''-' dw. {WA) 

Recognizing the integral as the Beta function B{k + j + N — k + I) = {k + j + j — 1)1 {N — k)\/ 
(7V + j + 7)!, 

Several comments are in order concerning {Dk{N)). First of all, although we restrict ourselves 
to discussing the first moment of Dk{N), we could in fact consider any moment {D^{N)) by 
taking [A~-^(ii))]" instead of A~^{w) in (IV.3). Doing so would alter 7 and the Cj's, but would 
not change our results qualitatively. Second of all, there is no loss of generality in taking our 
total surface area to be unity; scaling this area by a constant (or even, as might be more intuitive 
to statistical physicists, by N) would provide only a trivial scaling factor in our results. Third 
of all, we could imagine that the point x we consider is itself an {N + l)th site. This is simply 
a question of nomenclature. The problem of finding the expected distance from a point to its 
/cth-nearest site, for a system of N sites, is therefore equivalent to the problem of finding the 
expected distance hefween A'th-noarcst neighbors, for a system of N + 1 sites. 
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Let us now define the reduced variable {Dk{N)) by dividing out the leading asymptotic (large 
N) behavior from {Dk{N)): 

^ CO {k + j-iy. {N+j + j)V ^ • ' 

so that limjv^cx)(-C/!;-(A^)) = 1 (this is seen from StirUng's law). {Dk{N)) then provides the finite 
size scaling behavior for A;th-nearest neighbor distances. In what follows, we shall consider the 
properties of A~^{w), and its consequences on {Dk{N)). 



- 3 - Flat Surfaces 

On a flat surface, if we could neglect edge effects, the area included within distance I would 
simply be A(l) = ttP. In that case, A^^(w) = 1x1/1:, and so from (IV.3) and (IV.6) we would 
have 



The finite size scaling would thus be completely independent of k. 

As we are working with a surface of fixed (unit) area, however, we cannot avoid considering 
edge effects. Let us restrict ourselves to the case where the surface is everywhere locally Eu- 
clidean within some minimum neighborhood of radius Iq > 0. (The simplest example of this is 
a unit square with periodic boundary conditions, for which = 1/2- Obviously, many other 
constructions are possible.) Any required modifications to the A~^{w) expression in (IV.3) then 
concern only w greater than wq = A{lo). Correspondingly, (IV.4) remains valid up to remainder 
terms from the region of integration wq < w < I. Since the (1 — w)^~^''^ term in the integral 
is bounded above by (1 — wo)^~^~^ within this region, these remainder terms are exponentially 
small in N. Equation (IV. 7) is therefore still correct to all orders in 1/N, and may be written 
as the expansion 

(Dk(N)) = l- — + o(^ 

where all orders in are independent of k. The finite size scaling law for /cth- nearest neighbor 
distances on a 2-D flat surface without a boundary thus exhibits the universality in k to all orders 
in 1/N. 

The same holds true for flat manifolds in any dimension d. We assume there is an Iq such 
that the volume included within distance / < Zq is simply the volume of a d-dimensional ball: 

^^^^ = w °' 



A-^{w) = 



w 



-1 1/d 
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As before, the boundary conditions allow us to write {Dk{N)) up to remainder terms that are 
exponentially small in N, so from (IV.6), 

^ i/d + i/d' ^^f 1 



2N 

Thus for fiat spaces without a boundary, of any dimension d, the universality in A; holds to all 
orders in 1/N. 

It may be interesting to consider a slight variation on the problem, giving this universality 
exactly and not only to all orders. Take the case of a spherical surface embedded in 3-D Euclidean 
space, with the usual measure of area over the sphere, but with a peculiar sort of "distance" : 
rather than the conventional choice of the arc length (geodesic) metric, use the chord length. 
For a chord of length I originating at a pole of the sphere, the area of the spherical cap spanned 
by it is simply A{1) = nP. The /sth-nearest neighbor distance properties using chord length 
"distance" on this curved surface then appear analogous to those on a flat surface. There is, 
however, one important distinction. The relevant threshold wq for edge effects is in this case 
wq = 7r(2i2)^, where R is the radius of the sphere. Since 7r(2i?)^ is exactly equal to the total 
surface area of the sphere, it is set to 1. Equation (IV. 4) thus requires no corrections at all, and 
so the universality in (IV.7) is exact. 



— 4 — Curved Surfaces 

Let us turn to the case of a surface with true curvature, where the distance is now defined in 
terms of a metric, i.e., along geodesies of the surface. Let us consider again a spherical surface. 
With I now representing arc length, the area of the spherical cap spanned by an arc originating 
at a pole of the sphere is given by 



e 



1 — cos 



2R 



If the total surface area, A-kR"^, is normalized to 1, 



^(Osph e 



sin^ ■\/Trl, so 



A ^{■W)sphere = ^ sin ^ ^/w 

iw w (2j)! 



,7=0 



2j + 1 22i(j!)2 



(IV.8) 
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As in the case of the chord length "distance", this A ^{w) expression is exact everywhere for 
< w < 1. Equations (IV.3) and (IV. 6) then require no corrections, and we find 



j=0 



(2j)! (k+j -l/2)\ 



NNl 



2j + l 22i(j!)2 (A; -1/2)! (iV + j + 1/2)! 



1 + 



Ak-7 
24N 



(IV.9) 



Clearly, the k universality does not apply here. 

Another sort of universality, however, is found when wc turn to the more general case of 
an arbitrary closed surface, i.e., an abstract 2-D manifold with non-constant curvature and 
no boundary. For any surface, we may introduce a system of curvilinear coordinates u and v 
enabling us to write (at least piecewise) the differential length element ds in the conformal, 
orthogonal form (ElSENHART, 1909): 



ds"^ = f{u,v) [du^ + dv% 



(IV. 10) 



The Gaussian curvature K{u, v) of the surface may then be expressed in terms of the function 
f{u,v) by 



K 



2/3 



Of 



du 



Of 



dv 



it.) + if -f^-f 



(IV.ll) 



What is A{1) on this surface? In order to know this, we must first find out what are the manifold's 
geodesic lines. For ds given by (IV. 10), we may use the geodesic equation: 



ds2 ^ 2/ du 



[ds^ 



1 df du dv 

H = 0. 

/ dv ds ds 



(IV. 12) 



Let us expand u and v as functions of distance s from an initial point, along a fixed geodesic: 

o2 



U{s) = Uo + SUq + —Uq + 
V{s) = Vo + SVq + —Vq + ■ 



and 



(IV. 13) 



where uq = n(0), u'q = u'{0), etc., and likewise for v. Then, expanding f(u, v) in terms of u and 
V and substituting (IV. 13), 



/(u, v) = f{uo, Vo) + s [no/u(no, vo) + v^fuiuo, vo)] + 



u" v" 



I \2 



-fuuiuQ, Vo) + Uo^ofuviuQ, Vq) H T^fvv{uO, Vo) 



+ 0{s% 



where subscripts on / denote partial derivatives. 
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Using (IV. 10) and (IV. 12), we can solve for all but three of the coefficients in (IV. 13). Let 
us choose uq, vq and Uq to be these three. Now, consider the area A{1) about the point («o, vo). 
For ds given in (IV. 10), the differential surface element will be du = f dudv, so: 



Ail) 



f du dv 



I 



du dv 



ds du'o 

du dv 



ds du'n 



duQ ds 



ds du'n 



(IV. 14) 



The limits of integration over s are and Z; the limits of integration over u'q, which may be 
found from (IV. 10), are —^1/fiuQ, vq) and ^l//(uo, vq). In evaluating the Jacobian some care 
must be taken, as a sign ambiguity allows two solutions for the coefficients in (IV. 13). Ail) will 
be the sum of (IV. 14) evaluated at each of the two solutions, ultimately causing all odd powers 
of I to vanish. 

The result, after lengthy algebraic manipulations, may be written as 



Ail) = 



1 - 



ill 

24 p 



Ul + fv-fUu-fM + oii") 



(IV. 15) 



where in order to avoid cluttering the notation, we have omitted the (mq, vq) arguments at which 
all functions are to be evaluated. For the leading correction term in Ail) given by (IV.15), 
we recognize the expression (IV. 11) for the Gaussian curvature K. In retrospect, this is not 
surprising: by symmetry, the correction series to vrZ^ can contain only even powers of I, and if 
wc consider Ail) as a geometric expansion about a flat space approximation, K will be the only 
scalar quantity with dimensions of l^"^. 

With some perseverance, one may carry (IV.15) to higher orders in Z, and the coefficients of 
the expansion turn out to be 



Ail) 



nl' 



1 



iSK^ - 3[10(V/^:)^ + UKV'^K - hV^K]) + Oif) 



(IV. 16) 



161280 

where V is the gradient operator. We thus have a series expansion giving the area of a disc on 
any 2-D surface, in a form that depends only on intrinsic quantities, i.e., not on the choice of 
coordinate system. 

We may invert (IV. 16) to obtain the power series 



A-\w) 



K 9K'^ + 4V'^K 
^24^"^^ 19207r2 
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215047r3 



+Oiw^)] . 



(IV.17) 



Take as an example the special case of a spherical surface, where the Gaussian curvature is a 
constant K = l/R"^, or K = An foi a unit surface. All derivatives of K then vanish, leaving 



A (?i;)sphcro = Y r 



+ TTTT + Oiw^) 
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from which we recover the first few terms of our earlier result (IV. 8). 

Given an expression for A''^{w) on a general 2-D surface, we may find (i)fc(iV)) using (IV.3) 
and (IV.6). Let us consider the average J {Dk{N)) over the entire surface, which is found 
directly from the average of the series coefficients in (IV. 17). Now examine the 0{w) term in the 
series. By the Gauss-Bonnet theorem (ElSENHART, 1909), j K dji = 2ttx on any closed surface, 
where x is the Euler characteristic of the surface, a topological invariant. We therefore obtain, 
up to leading corrections, 



J A-^{w)dn 
j {Dk{N))dfi 



givmg 



NN\ 



l^X fc + 1/2 

(iV + 1/2)! [ 12 N + 3/2 

^ 2AN ^ V iV2 



(IV. 18) 



We thus discover another, very different sort of universality from the one we had in the case 
of fiat space. To 0{1/N), the finite size scaling law for fcth-nearest neighbor distances depends 
only on the surface's topology, and not on its detailed properties. 

The Euler characteristic x for a surface is related to its genus g hy x = 2(1 — g). Taking the 
torus as one example, y = 1, so x = and the ^-dependence in (IV. 18) once again disappears, 
at least to 0{1/N). This is to be expected: a flat space with periodic boundary conditions has, 
after all, the topology of a torus. And conversely, because of the topological invariance, all tori 
behave like flat space to 0{1/N). Taking the spherical surface as another example, g = 0, so 
X = 2 and we recover from (IV. 18) the result (IV. 9). 

The properties of J {Dk{N)) dfi are far less clear at higher order in 1/N. Using the divergence 
theorem and integration by parts, we may obtain from (IV. 17) 



J A-\w)dn = J 



K 



^ ^ 247r ^ 6407r2 



215047r3 



dn. (IV. 19) 



If we looked only at terms up through 0{w'^), we might believe that this scries is simply, by 
analogy with (IV. 8), an expansion of f{2/\/T{)sm^^ a/ Kw/Att d^. Unfortunately, starting at 
0{w^) we see this is not true, since the contributions of curvature and its gradients do not 
all vanish in the average over the surface! Furthermore, even for terms in (IV. 19) of the form 
J" d/x, at n > 1 there is no clear equivalent to the Gauss-Bonnet theorem; the theorem is a 
direct consequence of the integrand's linearity. Thus for a general 2-D surface, there does not 
appear to be a simplified form for the terms in J {Dk{N)) dfi beyond 0{l/N). More particularly, 
the only case in which J {Dk{N)) dji would be independent of k beyond 0(1/N) is if the curvature 
is identically equal to 0, i.e., a flat surface. 

Let us briefly consider the case of curved higher-dimensional manifolds. The calculation is 
now far more complicated, as it is no longer possible to write the metric tensor in a conformal 
form as we did in (IV. 10). In addition, whereas in 2-D the only intrinsic scalar quantity describing 
curvature is the Gaussian curvature K, for d > 2 there are d{d—l){d—2){d+3)/12 different such 
quantities (WeinberCx. 1972). However, all of them except K itself have at least dimensions of 
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It thus seems reasonable to conjecture that, as we argued in 2-D, the 0{P) correction term 
in A{1) can only involve K. (Indeed, we have verified that this is true in 3-D.) In that case, 
we may rely on the example of the spherical surface — easily generalized to d dimensions — to 
provide us with the initial terms for a general manifold: 



A{1) 
A-\w) 



.d/2 



{d/2)\ 
(d/2)! 
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d/2 



1/d 
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1/d 



(IV. 20) 



Note that A ^{w) now contains a series in w'^^'^ rather than in w. Appropriately modifying 
(IV. 3), it may then be shown that (-Dfc(iV)) is in general given by a series in for odd d, 

and l/Ar2/<i for even d. 

Consider, finally, the average J A~^{w)diJ. over the manifold. The higher dimensional gen- 
eralization of the Gauss-Bonnet theorem (Nakahara, 1990) involves an integrand of 0{l/l'^), 
or 0{l/w). The leading correction term from (IV. 20), J K dfi, therefore cannot be simplified 
further for d > 2; the only term that could possibly give rise to a topological invariant is the 
coefficient at 0{l'^), or 0{w). If d is odd, it is rather certain that no topological invariant will 
exist in the series. If d is even, the 0{w) term will first contribute to the J {Dk{N)) dji series at 
0{1/N) — as in 2-D, although at higher dimensions this will no longer be the leading correction 
term. While we cannot rule out the possibility of indeed obtaining a topological invariant at 
0{l/N), the 0{w) term in A-^{w) is in general a complicated one involving many different 
curvature scalars, so this is far from obvious. It remains an open question. 



— 5 — Regge Calculus 

We have remarked that from a physical point of view it is natural, in the 2-D case, for the 
leading corrections in A{1) to contain only the Gaussian curvature, as this represents the leading 
deviation from planarity. Consequently, only the mean curvature — or, using the Gauss-Bonnet 
theorem, the Eulcr characteristic x — matters in the 0{1/N) term of /" {Dk{N)) d^. We have 
seen using differential methods (geodesies) that this physical picture is indeed correct. These 
methods apply to a smooth surface. For polyhedral surfaces, which are not smooth, we may 
in fact obtain a similar result more easily, using the non- differential method of Regge calculus. 
Consider a polyhedron with a number of vertices, edges, and faces. Following the work of Regge 
(Regge, 1961) and others since then (Cheeger, Muller & SCHRADER, 1984), we observe 
that the curvature is concentrated at the vertices and is measured by a deficit angle: if 9i is 
the sum of the angles incident on vertex i, the deficit angle at that vertex is Aj = 2n — 9i. It 
may then be shown that the Gauss-Bonnet theorem, on polyhedra, reduces to Euler's relation 
27rx = EiAi. 

Let P be a polyhedron with a fixed number of vertices, and consider the problem of finding 
the finite size scaling J {Dk{N)) djJL on P. As N 00, corrections to the flat space value about 
a given point x arise only when x is near one of the vortices, because only then can curvature 
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(i.e., tiie deficit angle) enter into the local calculation of A{1) about x. It is then sufficient to 
understand the corrections associated with one vertex at a time. Consider a vertex i. A{1) 
receives a contribution from i that depends exclusively on the deficit angle Aj. For small deficit 
angles (small corrections to A{1)), a linear approximation may be used and this contribution 
will simply be proportional to Aj. Correspondingly, the leading corrections to A{1) — and thus 
the 0{l/N) term in the finite size scaling series — will be proportional to Aj. Summing over 
all the vertices i = 1, . . . , A^, we find that the leading correction term in J {D^^N)) dfi is indeed 
a linear function of %, and we recover the topological invariant derived in the case of a smooth 
manifold. 

A word of caution is necessary, however. It is tempting at this point to take the limit where 

P becomes a smooth manifold, thus recovering (IV. 18). Unfortunately this will not work; a 
direct computation shows that the limit does not commute with the limit N ^ oo taken above, 
and the coefficient thus obtained for the 0(1/N) term will not be the correct one. 

— 6 — Conclusions 

We have considered the finite size scaling of mean distances to neighboring sites, when N sites 
are distributed randomly and uniformly on a surface with no boundaries. When the surface 
is flat, we have found that the entire 1/N series describing the mean fcth- nearest site distance 
is independent of k. This universality applies equally well to higher moments of the distances, 
and to Euclidean manifolds in dimensions greater than 2. For surfaces with curvature, while 
this general property is no longer valid, we have found that the leading correction term in the 
series, averaged over the surface, is a topological invariant. The scaling series thus depends, to 
0{1/N), on the genus of the manifold but not on its other properties. 

Although we have considered these universalities only for the moments of point-to-point 
distances, similar properties hold for higher order simplices such as areas of triangles associated 
with nearby points. The problem is thus a natural one to consider further in the context of 
random triangulations, foams or other physical problems (Itzykson & Drouffe, 1989) tightly 
connected to geometry. 
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Computational complexity: P vs. NP 

The inherent challenge of the TSP lies in the fact that no known algorithm can find the optimal 
tour of an arbitrary A^-city instance in a number of steps polynomial in N. In computational 
complexity theory, this is what classifies it as a "hard" problem. An "easy" problem would be 
one that could be solved by a polynomial algorithm. The word "easy" , of course, does not imply 
"fast": an 0{N^) algorithm, at N = 10,000, could already involve immense running times! 
There is nevertheless a large conceptual difference between a problem that can be solved in no 
worse than polynomial time, and a problem that requires, say, exponential time. 

Computer scientists formalize these concepts with the notion of P and NP. A problem belongs 
to class P if there exists an algorithm that solves the problem in a time growing polynomially, 
or slower, with the size N of the problem. A problem belongs to class np if it is merely possible 
to test, in polynomial time, whether a certain "guessed" solution is indeed correct. Note that 
NP stands for "non-deterministic polynomial", and not, as is sometimes mistakenly thought, for 
"non-polynomial". P is in fact a subset of np, since for any problem whose solution can be found 
in polynomial time, one can surely verify the validity of a potential sohition in polynomial time. 

Formally speaking, the classes P and NP apply only to decision problems, where the solution is 
simply a "yes" or a "no" . It is however, generally possible to relate combinatorial optimization 
problems to associated decision problems. For the case of the TSP, instead of asking for the 
optimal tour length, we may phrase the associated TSP-decision as follows: is there a tour whose 
length is less than a given value B? It is relatively straightforward to see that the optimal tsp 
tour length can be found to arbitrary precision by executing a sequence of TSP-decisions (with 
bound B, say, decreasing in discrete steps), where the length of the sequence is polynomial in N 
(Johnson & Papadimitriou, 1985). The TSP-decision is thus said to be polynomially reducible 
to the TSP; this implies that the two problems, though not necessarily themselves equivalent, 
are at least of equivalent computational complexity. 

In order to see that the TSP-decision is in NP, let us use the following more careful definition 
of the NP class (CoOK, 1971). A decision problem belongs to NP if any instance calling for a 
"yes" response contains a certificate, itself of size polynomial in N, allowing the "yes" answer to 
be verified in polynomial time. For the TSP-decision, this certificate is simply a tour of length 
less than B: the certificate is of size N, and in 0{N) steps one may verify that it is a legal tour 
and satisfies the length bound. The certificate is therefore the "guess" that is tested. Of course, 
we cannot hope to verify a "no" answer in the same way as a "yes" answer; the only property 
characterizing a "no" instance is that no certificate exists. The fact that a decision problem is 
in NP merely means that we can confirm a "yes" instance if we happen to be given the right 
certificate. 

In the case of the TSP-decision, one can go further. An important subset of the class of 
NP problems is the class of NP-complete problems. These are problems that are as least as 
complex as any other np problem — in other words, an algorithm capable of solving an NP- 
complete problem could be mapped onto an algorithm capable of solving any NP problem, via a 
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polynomial encoding. (This is not quite the same thing as the notion of polynomial reducibility, 
as the encoding associates a single instance of one problem with a single instance of another, 
rather than with a sequence of instances.) It has been proven that the TSP-decision is NP- 
complete (Papadimitriou, 1977)- Let us consider the implications of this. We have already 
noted that P C np, since a problem solvable in polynomial time has by definition a certificate 
(its solution!) verifiable in polynomial time. If it turns out that a polynomial algorithm could 
be found for solving the TSP-decision — or any other NP-complete problem — one would exist 
for all NP problems and we would have P = NP. No such polynomial algorithm has ever been 
found, and it is conjectured that none exists, so that P 7^ NP. Proving this conjecture, however, 
has been an open question in complexity theory since the 1970s. 

Note that although it has become standard practice to refer to the TSP as being NP-complete, 
the more correct term is actually NP-Ziard. This is because, strictly speaking, the complexity 
classes P, np and NP-complete are reserved for decision problems. NP-hard is the more general 
designation for problems to which NP-complete decision problems are polynomially reducible. 
The TSP-decision is NP-complete; the TSP itself, being of equivalent complexity, is thus NP-hard 
(see Garey & Johnson, 1979; Johnson & Papadimitriou, 1985). 

A different issue from that of solving an NP-hard problem in polynomial time is that of 
approximating its solution in polynomial time. For the Euclidean tsp, Karp (1977) developed 
the fixed dissection algorithm, which generates a tour using the "divide-and-conquer" strategy of 
partitioning space into subparts, solving the TSP for the cities within each subpart, and joining 
up these subtours in such a way as to form one large (non-optimal) tour through all cities. Karp 
proved that this algorithm can with high probability (prob. ^ 1 as ^ 00) give tours whose 
length is within (1 -|- e) times optimality, for arbitrary e. The expected execution time of this 
algorithm is 0{N'^ logN). 

Karp's construction has three limitations: first, the heuristic does not guarantee the 1 + e 
bound at finite N; second, the running time could in a worst-case situation be exponential; third, 
the method only applies to instances from the random ensemble, with a uniform distribution 
of cities. Recent work by Arora (1997) appears to have resolved these three difficulties in 
the case of the Euclidean tsp. Using a different partitioning scheme, Arora's algorithm finds a 
guaranteed (1 + e)-approximation in 0((log A^)'^^'^/^^'' ^iV) time, for arbitrary instances in d 
dimensions. This sort of method, unfortunately, cannot easily be generalized to other problems. 
It has in fact been proven (Arora, Lund, MoTWANi, Sudan & Szegedy, 1992) that if p 7^ 
NP, (1 -|- e)-approximation algorithms cannot exist for all NP-hard problems. The fact that one 
has been found for the Euclidean TSP is quite remarkable, and suggests that further refinement 
of the NP-hard classification may be necessary. 
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Outline of a self-averaging proof 



The original proof of self-averaging in the Euclidean optimal tour length LE{N,d), by Beard- 
WOOD, Halton & Hammersley (1959), is quite technical. There is however a more accessible 
proof, by Karp & Steele (1985), using the fixed dissection algorithm (Karp, 1977) mentioned 
in Appendix A. This algorithm generates a (non-optimal) tour by dividing the space into sub- 
parts, finding the optimal subtour within each subpart, and connecting up all of these subtours 
(minus one link, in each case) to form one large tour. Karp and Steele also took advantage of 
a lemma (proved by induction), stating that in the d-dimensional unit hypercube there must 
always exist a tour of length less than dN^~^^^ + ddN^~^^^^~^^ {Sd is a constant depending only 
on d). 

They then considered a Poisson process placing points with unit intensity in the hypercube 
[0, t]'^, and looked at the expectation value F{t) of the optimal tour length through these points. 
If {LE{N,d)) is the expected optimal tour length through N points in the unit hypercube, t 
times this quantity will be the expected optimal tour length through N points in [O,*]''. Prom 
the Poisson distribution, we thus obtain: 

m=T.^~' -j^t{LE{N,d)). (B.l) 

N=0 

Now, the fixed dissection algorithm bounds the optimal tour length through the cities placed 
by the Poisson process. This upper bound is equal to the lengths of the various subtours plus 
the length of the large circuit through space needed to connect the subtours. Bounding the 
latter length using the lemma above, and partitioning space into m'^ equal subparts, we then 
obtain the expectation value bound: 

F{t) < m'^F(t/m)+i(d(m'^)^-i/'^ + 5d(m'^)^-^/('^-^)), or (B.2) 
m ^ F{t/m) , d m-^jd-^^ 

Clearly F{t) is monotone increasing, and because of the bounding lemma, F{t)/t'^ must be 
bounded. As t — 00, then, for fixed m, (B.3) implies that F{t)/t'^ approaches a limiting value, 
which we shall call PE{d)- Prom (B.l), substituting u = t'^, 

^ {LEiN,d)) _ 

N=0 

From the Tauberian theorem (see Karp & Steele, 1985, p. 187) for the Poisson distribution, 
finally, this gives: 
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The statement is in fact stronger than it appears: Karp and Steele are also able to 
bound the variance of LE{N,d) using an inequality due to Efron & Stein (1981), giving 

now know, therefore, not only that the distribution of 
LE{N,d)/N^~^/'^ becomes increasingly sharply peaked at large N, but also that its width goes 
to zero as l/\//V, i-e., as a Gaussian. 
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Estimates for a non-uniform distribution of cities in space 



We have seen that the optimal tour length is self-averaging, in the ensemble of cities indepen- 
dently and uniformly distributed in Euclidean space. With probability 1, the random variable 
LE{N,d)/N^~^/'^ tends to a constant l3E{d) in the N ^ oo limit. The original proof of this 
result, due to Beardwood, Halton & Hammersley (1959), is however more general and 
deals with an arbitrary density p(x) dx. of independently distributed cities in space. The general 
statement of self-averaging is that, with probability 1, 

where f3E{d) is independent of the density, and thus equal to its value for a uniform distribution. 
What this means in practice is that the numerical value of Psid) is of relevance for cities 
distributed independently in space with any distribution. 

In this appendix we propose to use this result to give an a priori estimate oi Le for any 
instance. Note, first, that the factor J[p(x)]^~^/'^ dx applies just as well to quantities such as 
the distance Di{N,d) between nearest neighbors. For an arbitrary Euclidean instance at large 
N, then, we may estimate the optimal tour length ^non-umform £qj. ^-^^^ instance as: 

/ r uniform \ 

^non-uniform ^ ^non-uniform \ E /_ 

E ~ 1 ^^uniform^ ^ ' 

where DJ^on-uniform ^Yie mean nearest neighbor distance for the instance we are considering. 
(L^"'*^"^™) is given by its value in the Euclidean uniform ensemble, and ^/)u°if°r™^ ig the ensemble 
average of nearest neighbor distances. Both of these quantities, as functions of N and d, have 
been calculated in Chapter II. Note that there is no need to use units where volume equals 1 
in the non-uniform instance under consideration; our use of D^on-uniform 2) will take into 

account any volume scaling. 

Let us try this out on a tsp instance well known in the operations research literature, 
the AT&T-532 instance. -"^ This consists of 532 sites belonging to the U.S. telecommunications 
company. For clear commercial reasons, the distribution of sites in the instance follows roughly 
the population density of the country, and thus is highly non-uniform. The instance was first 
solved to optimality (see Figure C-1) by Padberg 8z RiNALDi (1987), and in their units, the 
optimum was found to be 27,868. How closely does our approach predict this result? In their 
units, once again, the mean distance between a city and its nearest neighbor can be measured 
as 36.205. Using the expressions given in Chapter II, we would then estimate: 



^AT&T-532 is one of the many instances found in TSPLIB, a library assembled by Reinelt (1991) and available 
from http://softlib.rice. edu/softlib/tsplib. 
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/ r uniform \ 

^non-uniform ^ ^non-uniform \ E /_ 

E ~ 1 ^^uniform^ ^ ' ' 

« £)non-uniform ^ ^^^3) X (^1 - - X 2iV (C.4) 

^ 27,427 (C.5) 

which is within 1% of the actual optimum. 

It may seem surprising that our finite size scahng law given periodic boundary conditions pro- 
vides such a good estimate for a tour length calculated given open boundary conditions. There 
are two reasons for this agreement. First of all, boundary effects are involved in £)non-uniform 
(C.2) incorporates these effects into the i^non-umform estimate, and although they may not be 
precisely the correct effects for the optimum tour length, they are undoubtedly close. Second of 
all, (3e{2) itself docs not depend on the boundary conditions chosen. For a discussion of why this 
is so, the reader is referred to JAILLET (1993), who has proven that open and periodic boundary 
conditions give the same Psid) value. 
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Numerical methodology 

In this appendix we discuss the algorithms used to obtain our simulation results. For both the 
Euclidean and random link TSPs, we have performed runs at instance sizes N = 12, 13, 14, 15, 
16, 17 using the LiN & Kernighan (1973) heuristic (lk), and at = 30 and N = 100 using 
the Chained Local Optimization heuristic (CLO) of Martin, Otto & Felten (1991). 

LK heuristic 

The LK algorithm may be sketched as follows. The kernel of the algorithm is what might be 
termed LK-search. We start off with an initial (non-optimal) tour. LK-search takes an arbitrary 
starting city — for convenience, we generally start with the top city on a fixed "list" in our 
program. Call this city iQ. Now pick a starting direction — either "forwards" or "backwards" 
— and call ii the next city in that direction along the initial tour. Call li the link between 
io and ii, and remove that link. We will now attempt to reconnect ii to a new city i[ that is 
somewhere else on the tour, resulting in the situation shown schematically in Figure D-1. Let 
I'l be the link to that new city i'l- How do we choose i\, and thus l[? i\ is the nearest neighbor 
to ii that was not already connected to ii in the initial tour. There is another very important 
requirement for i'^. it must be such that ![ < li. This it the gain criterion, and is applied at 
all stages of LK-search. If the gain criterion cannot be met, the search using this io and is 
abandoned. If the gain criterion is met, however, the search deepens. 




Figure D— 1: One step of LK-search, showing the removed link li (dashed line) and added link l[ (bold 
line). 
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Figure D-2: LK-search at step 2, showing tadpole with two hnks removed and two Unks added. 



At this point we have what might be called a "tadpole graph" (the more technical term is 
a one-tree). We now pick one of the two cities next to i'^ along the initial tour, call it Z2, and 
remove the link I2 joining them. Which of the two possibilities do we use for 12, and thus ^2? We 
pick ^2 such that, if ^2 were then to be connected to the dangling end (io), a single closed tour 
would result. Now continue as before. We attempt to reconnect ^2 to a new city i'2, resulting in 
the situation in Figure D-2. I2 is the link to the new city, and that city, as before, is the nearest 
neighbor that was not an adjoining member of the original tour. The gain criteria get stronger 
as the search gets deeper: the total gain must remain positive, so we must have I'l + I2 < h +I2, 
or else the search aborts. LK-scarch continues recursively in this way, with the vertex of the 
"tadpole" hopping around with the end point io staying fixed, until the gain criteria force it to 
abort. There is, however, an important additional part to the gain criteria. At each step m, 
whilst im is being chosen in such a way as to make it possible to close up the tour, the length 
of that closed tour is recorded. If at any step of LK-search the tadpole length for that step is 
longer than the best closed tour recorded so far, the search terminates. 

When a search terminates, one of two things can happen. If the best recorded tour is shorter 
than the initial tour, the entire process begins again but using that tour as the new onc.^ If, on 
the other hand, the best recorded tour is not shorter than the initial tour, we recommence with 
the alternate choice of ii: if the starting direction was "forwards", now we try "backwards", 
and vice versa. When, however, both choices of ii have been exhausted, we select a new io by 
advancing to the next city on the list. The algorithm, in our implementation, terminates when 

^This will happen, notably, if the search is abandoned because of the second part of the gain criteria (closed-up 
tour shorter than tadpole), as the first part of the gain criteria (J^^ I'i < k) requires that a tadpole always be 

shortcT than tin' initial tour. 
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Figure D-3: Double-bridge change executed by CLO heuristic, generated by removing links shown with 
dashed lines, and reconnecting them differently. 

an entire "pass" of N starting cities io results in no improvement. 
CLO heuristic 

CLO is a stochastic algorithm, combining large Monte Carlo jumps with embedded LK local 
search. In our implementation, it works as follows. A (non-optimal) initial tour is optimized 
using LK. This brings it to an LK local minimum. An attempt is then made to modify the tour by 
a random 4-change (4 bonds disconnected and then reconnected differently), which would not be 
accessible by the kind of sequential changes performed by LK-search. The 4-change in question 
is known as a double-bridge kick, shown in Figure D-3: one 2-change disconnects the tour, and 
the other 2-change reconnects the two parts in a different place. After this 4-change is carried 
out, LK is again used to optimize the resulting tour. If the new local optimum is better than 
the previous one, the attempt succeeds, and CLO iterates the random 4-change procedure from 
the new tour. If not, the attempts fails and CLO tries another random 4-change on the old tour. 
This continues for a fixed number of steps. The idea is that, by directing its search through state 
space, CLO performs better than simply running LK from random starts an equivalent number 
of times. 

Use of heuristics 

For instance sizes in the range 12 < iV < 17, we used the LK heuristic. Our method consisted 
of two parts: first, we used a testbed of instances to estimate the systematic bias in lk arising 
from the fact that LK does not always find the true optimum; second, we did our "production 

runs" to determine the actual optimal tour lengths. 



93 



Appendix D — Numerical methodology 



In order to determine the systematic bias, we performed the following procedure on values of 
N from 12 through 17: (i) we generated a testbed of 200 random instances; (ii) for each random 
instance, we generated 100 different random starting tours; (iii) for each random starting tour we 
ran the LK heuristic. Based on the assumption that the best LK-opt obtained over 100 random 
starting tours was indeed the true optimum for that instance, we calculated the expected bias 
per instance. We then averaged this bias over the testbed of instances, obtaining an estimate 
for each value of A'^. 

For our production runs, we then performed the following procedure on values of N from 
12 through 17: (i) we generated 100,000 random instances; (ii) for each random instance, we 

generated 10 different random starting tours; (iii) for each random starting tour we ran the LK 
heuristic. Wc then took, for each instance, the best LK-opt obtained over the 10 random starting 
tours, and averaged over all instances, obtaining an estimate of the mean optimum length for 
each value of N. 

For instance sizes iV = 30 and N = 100, we used the CLO heuristic with 10 steps. The method 
and procedures were identical to those for LK, apart from the numbers for the production runs: 

for = 30 we generated 8,000 random instances and 5 random starting tours per instance; for 
N = 100 we generated 1,200 random instances and 20 random starting tours per instance. 

Other operational issues 

The choice of number of random starting tours per instance in our production runs was motivated 
by the need to keep the systematic bias to a minimum. For the LK runs, using 10 random starting 
tours per instance allowed us to have in the worst case, at N = 17, a bias of 1 part in 200,000 
for the Euclidean tsp and 1 part in 170,000 for the random link tsp. This bias was therefore 
negligible at the level of numerical precision used (4 decimal places). For the CLO runs at 
N = 30, using 5 random starting tours per instance was sufficient to give us an estimated bias 
of under 1 part in 400,000 for the Euclidean tsp, and under 1 part in 800,000 for the random 
link TSP. We suspect these figures — particularly the latter — might be too good to be true, 
so it is possible that our estimate is corrupted somewhat by instances in the testbed where the 
true optimum was never found. (Even if our estimate of the bias is off by an order of magnitude, 
though, it will still be negligible compared with our statistical error bar, which at = 30 is 
about 1 in 3,000 for both Euclidean and random link.) For the CLO runs at = 100, using 20 
random starting tours per instance gave an estimated bias of 1 part in 120,000 for Euclidean 
and 1 part in 60,000 for random link — again, comfortably below the statistical error of 1 in 
6,000 for Euclidean and 1 in 2,000 for random link. 

Finally, let us note a few further details concerning our code. Timings for our d = 2 random 
link simulations are given in Table D-1; as can be seen we ran on four different machines, two 
of them Dec Alphas and two of them Sun sparcs. The code itself is written in C, and is a 
heavily modified version of the CLO package developed by Steve Otto (otto@cse.ogi.edu) and 
Robert Prouty (prouty@cse.ogi.edu). In the code, cities are placed on a 10,000 x 10,000 grid, 
and distances are rounded to the nearest integer. We confirmed that the mesh was fine enough 
for roundoff error to be negligible to within the 4-decimal-place precision of our numerical data. 
The random number generator Tiscd. from the CLO package, is a standard linear congxTiential 
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Table D-1: Simulation timings for c? = 2 random link. Instance sizes N = 12 through N = 17 used LK, 
N = 30 and N = 100 used 10-step CLO. # runs indicates number of instances times number of random 
starts per instance. 



N 


lk/clo 


# runs 


Machine 


CPU time (sees) 


12 


LK 


500,000 


Dec Alpha200 


7133.4 


12 


LK 


500,000 


Sun Ultral 


22001.2 


13 


LK 


500,000 


Dec Alpha200 


10772.0 


13 


LK 


500,000 


Sun Ultral 


26025.2 


14 


LK 


500,000 


Sun SparclOO 


32755.3 


14 


LK 


500,000 


Sun SparclOO 


32420.9 


15 


LK 


500,000 


Dec Alpha200 


10724.9 


15 


LK 


500,000 


Dec AS 1000 


9357.5 


16 


LK 


500,000 


Sun SparclOO 


41459.0 


16 


LK 


500,000 


Sun SparclOO 


38795.3 


17 


LK 


500,000 


Sun Ultral 


43508.3 


17 


LK 


500,000 


Sun SparclOO 


46584.2 


30 


CLO 


5,000 


Sun Ultral 


10801.6 


30 


CLO 


5,000 


Sun SparclOO 


11717.5 


30 


CLO 


5,000 


Sun SparclOO 


11621.3 


30 


CLO 


5,000 


Sun Ultral 


11163.9 


30 


CLO 


5,000 


Sun SparclOO 


10773.3 


30 


CLO 


5,000 


Sun Ultral 


10569.0 


30 


CLO 


5,000 


Sun SparclOO 


11558.9 


30 


CLO 


5,000 


Sun SparclOO 


11287.7 


100 


CLO 


6,000 


Dec Alpha200 


19694.6 


100 


CLO 


6,000 


Dec AS 1000 


17982.2 


100 


CLO 


6,000 


Sun SparclOO 


71671.1 


100 


CLO 


6,000 


Sun SparclOO 


71137.9 
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algorithm working as follows: 

#dcfinc MASK ( OxTfffFfff ) 
#definc MULT 1103515245 
#define ADD 12345 
#define TW0T031 2147483648.0 
AAA = MULT & MASK; 
BBB = ADD & MASK; 

where & is the binary AND operator. The routine then updates the integer quantity randx to 
be (AAA*randx + BBB) & MASK, and returns a double-precision variable randx/TW0T031. 
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Another numerical study of /3{d) 

Subsequent to our work, Johnson, McGeoch & Rothberg (1996) have used a variant of our 
approach to confirm the values /3e{2) = 0.7120 ± 0.0002 and /3e(3) = 0.6979 ± 0.0002 (to higher 
precision) and to extend the results to d = 4. Let us briefly summarize their work, and the 
differences between our approach and theirs. 

Our finite-size scaling analysis, presented in Chapter II, uses data points at values of N from 
N = 12 to N = 17, as well as the points N = 30 and N = 100. Johnson et al, making use of 
our scaling law (II.2), work with larger instances sizes — N = 100, N = 316 and N = 1,000 — 
where a better fit may presumably be obtained. This requires much more powerful computational 
resources than our own, but also good algorithms and efficient methods for reducing statistical 
error to a minimum (so as to avoid too large a number of time-consuming runs). 

When one wishes to estimate {LE{N,d)), the most direct way is to take the numerical 
average Le{N, d) over a sample of K randomly chosen instances. This estimator has an expected 
statistical error (t{K) = ai^^l\fK^ where ctl^ is the instance-to-instancc standard deviation of 
Le- One can improve on this, however. Let us instead measure Le{N, d) — L*{N, d) + L*(N, d). 
If L* is a quantity closely correlated with Le, Le — L* will have a substantially lower variance 
than Le itself; if in addition L* has a small variance (or better yet, if we can measure the 
ensemble average (L*) exactly) then we will ultimately obtain a measurement of (Le) to higher 
accuracy than with the direct estimator Le- This is the method we have adopted in our own 
work. For L*, we use A(Li + L2)/2, where A is a parameter, Li is the mean distance between 
nearest neighboring points in an instance, and L2 is the mean distance between second-nearest 
neighboring points in an instance. The great advantage of this approach is that an analytical 
expression can be found giving (Li + L2)/2 exactly, so the only statistical error is due to the 
estimator Le — A(Li + L2)/2. A can be chosen so as to minimize the variance of this estimator: 
it is relatively simple to show that the value A* minimizing A is given by 

(Le(Li + L2)/2)-(L^)((Li+L2)/2) 
A = J- ^2 (b.lj 

where cr(^i^_^_i^y2 is simply the instance-to-instance standard deviation of (Li + L2)/2. Notice 
that A* is in fact a measure of how closely correlated Le and (Li -|- L2)/2 are — if they are 
perfectly correlated. A* will be equal to 1, and if they are uncorrelated. A* will be equal to 0. 
In practice, we have found that the degree to which they arc correlated varies little in A'^. For 
the d = 2 Euclidean case, we found A* ~ 0.75 at = 15; keeping A* fixed to this value over all 
A*", the variance reduction ranged from 0.38 at = 15 to 0.43 at N = 100. We thus succeeded 
in reducing the overall statistical error a{K) by a factor of between 2.33 and 2.63. (For the 
random link TSP, the same method reduced the error by a factor of between 2.69 and 3.08.) The 
relative insensitivity of the variance reduction scheme to N suggests that correlations between 
Le and (Li + L2)/2 are stable at large N; a related property has been seen in the random link 
case in Chapter III (see Figure III-5). 
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In the work of JOHNSON, McGeoch & Rothberg (1996), a somewhat different procedure is 
used for reduction of variance. The emphasis of their article is on analyzing the Held-Karp lower 
bound, the solution to an integer programming relaxation of the TSP (Held 8z Karp, 1970). 
Thus instead of measuring (Lii) using Le — L* + L*, they use the estimator L^/Lhk Lhk, 
where Lhk is the Hcld-Karp bound. The advantage of using this quantity is that it turns out 
to be correlated very closely with Le, so that Le/Lhk has extremely low instance-to- instance 
fluctuations. (To compare: at A'' = 100 in 2-D, our estimator Le — L* leads to a statistical 
error for sample size K of about 1.3%) /s/K, whereas theirs leads to a statistical error of about 
0.3%/-\/K.)^ Furthermore, in measuring Le/Lhk Johnson et al. used exact codes for the 
A'^ = 100 and N = 316 instances, so there is no systematic error in these measurements (for 
= 1,000 they used the same large N heuristic as we used for = 30 and A^ = 100). For the 
second part of the estimator, Lhk itself, they used a rapid approximate method — though they 
were able to quantify the systematic bias inherent in this approximation (about 0.005%) on the 
basis of a testbed of instances where Lhk had already been found exactly. They then corrected 
for this systematic bias. 

Note that while in our measurements, (L*) is calculated exactly over the ensemble and 
the statistical error is due exclusively to the fluctuations of L^; — L*, in the measurements of 
Johnson et al., Lhk involves a significant statistical error whereas Le/Lhk involves an almost 
negligible error. It is instructive to compare our A^ = 100 results with theirs (this is the only 
2-D instance size we have in common with them). Our data, on the basis of 6000 instances, 
give (Le(100,2) - L*(100,2)) = 2.4342 ± 0.0012 and (L*(100,2)) = 4.6934, so (Le(100,2)) = 
7.1276 ± 0.0012. Johnson et al. obtain {Le{100,2)/Lhk{W0,2)) = 1.005542 ± 0.000027 on the 
basis of 13,957 instances, and {Lhk{WO, 2)) = 7.0897 ± 0.0006 on the basis of 98,246 instances, 
giving (Le(100, 2)) = 7.1290±0.0008. We are unaware of the running times for the latter results, 
and thus cannot compare the relative efficiency of the two approaches. It is clear, however, that 
our resources could not permit instance sizes much in excess of A/" = 100, whereas their tests 
went up to AT = 1,000. 

There is, interestingly, a further bias in the estimator used by Johnson et al. that they have 
not noted. This arises from the fact that, of course, {Le) 7^ {Le/Lhk) {Lhk)- It is relatively 
straightforward to calculate this bias at large Le and Lhk- 

I \ ,j \^/r \ / 1 , '^Lhk {LeLhk) - {Le){Lhk) , \ „x 

\ w / ^"^^^^ ^ ^"^^^ I,' + (w {Teht^) — ■ ■ ■ J ^^-'^ 

where o"/,^^ is the instance-to-instance standard deviation of Lhk- In 2-D, the bias here is of 
0{1/N). Without taking account of this bias, Johnson et al. use our scaling law — truncated 
to subleading order — to fit data at N = 100, N = 316 and N = 1,000, obtaining /3b (2) = 
0.7124 ± 0.0001. Fortunately, the bias affects only the subleading terms {0{1/N) and beyond) 

^The estimator of Johnson et al. follows a similar method used earlier by Sourlas (1986), who measured 
LE/Lth Lth, where Lth is a weighted average of fcth-nearest neighbor distances for k up to 5. (This is by contrast 
with our own L*, which is simply an unweighted average of k up to 2.) Since Sourlas' results were for the d = 1 
random link case, it is difficult to compare with our simulations; let us note, however, that at iV = 100 he found 

his method reduecxl (t(K) by a factor of 4. 
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in the fit and should not affect the leading term that gives Pe (though care must be talcen when 
using their Le results for finite N). What they do not consider, however, is a test of goodness- 
of-fit; instead they estimate the I3e{2) error bar "conservatively" by adding up the error bars 
for the three data points in the fit.^ 

Let us attempt to perform the analysis somewhat more carefully, in order to provide a 
more meaningful comparison between their conclusions and ours. As L^; and Lhk clearly are 
correlated very closely, we shall make the assumption that the bias in (E.2) can in fact be 
neglected — at least &t N = 316 and N = 1,000 — and simply combine their data at these two 
points with our own data at smaller N. Table E-1 sumarizes the effects of this, comparing the 
coefficients in the scaling scries {Le^N, 2)) = ^/N /3e{2) [1 + A{2)/N + B{2)/N'^], without and 
with these two new data points. The two results are consistent with each other: the coefficients 
of the fit are relatively stable (even, surprising 

ly to 0(l/iV2))^ and the new is consistent with 
7 degrees of freedom (10 data points minus 3 fit parameters). The error in /3e(2) is obtained by 
the standard procedure of determining the values of this quantity that make increase by 1; 
those values are then /3£;(2) plus or minus one standard deviation. 

Table E-1: Values of coefficients of ffi {Le{N,2)) = VN (3e{2) [1 + A{2)/N + B{2)/N% for previous fit 
{N = {12,13,14,15,16,17,30,100} and new fit {N = {12,13,14,15,16,17,30,100,316,1000}). 



Result (3e{2) A{d) B{d) ^ 

Previous fit (8 data points) 0.7120 ± 0.0002 0.1088 -1.064 5^ 
New fit (10 data points) 0.7123 ± 0.0001 0.0982 -0.9819 7.97 



If we stretch our assumption further and assume that the bias in (E.2) is negligible even 
at = 100, and therefore combine our own N = 100 data with that of Johnson et al. — 
obtaining (Le(100,2)) = 7.1286 ± 0.0007 — the results of "New fit" do not change at aU, to 
the precision shown! Only changes, decreasing now to 7.03. Another possibility is to fit 
exclusively to their 3 data points, though one must be cautious when using results based on so 
little information. Nevertheless, a two-parameter fit of the form \/N I3e{2) [1 + A{2)/N] gives 
I3e{2) = 0.7123 ± 0.0001 and A{2) = 0.0791 (with = 0.016, though this value is purely 
anecdotal as there is only one degree of freedom!). Thus, this (3e{2) result indeed appears 
credible. It is unclear why it is quoted in their article as 0.7124; however, assuming that the 
values quoted for the corresponding results at d = 3 and d = 4 are correct. Table E-2 compares 
our results and theirs. 

Finally, Johnson et al. also performed simulations on the d = l random link case, although 
they did not attempt to extrapolate to /3^i(l), instead noting simply that their large N data were 
consistent with the asymptotic value 1.0208.^ Let us perform a fit of the usual sort on their data, 

^Furthermore, it appears that in the notation of Johnson et al., the error bar indicates the extremes of the 
95% confidence interval, hence ±2a. In quoting their results, we use the more standard notation of ±cr. 

^Their notation differs from ours by a scaling factor of 2; we cite their results following our own notation, 

without this factor. 
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Table E— 2: Comparison of Psid) from two different numerical studies. 

d (3E{d) from our study f3E{d) from Johnson's study 
~2 0.7120 ± 0.0002 0.7123 ± 0.0001 

3 0.6978 ± 0.0002 0.6980 ± 0.0002 

4 N/A 0.7234 ± 0.0002 



using the 5 points that they give: N = 100, N = 316, N = 1,000, N = 3,162, and N = 10,000. 
(They do not actually carry out simulations of Lri/Lhk aX N = 10,000, considering that at the 
level of precision used this quantity is indistinguishable from 1, based on its value of 1.000036 at 
= 3,162.) The results are plotted in Figure E-1: I3rl{1) = 1.0209 ±0.0002, and = 2-62 for 
2 degrees of freedom (5 data points minus 3 fit parameters). This provides excellent experimental 
confirmation of the cavity predictions at d = 1. 



1.03 




1.01 



0.0025 



0.005 
UN 



0.0075 



0.01 



Figure E-1: Finite size scaling oi d = \ optimum. Best fit (x^ = 2.62) is given by: {L]il{N,1)) = 
1.0209(1 - 0.1437/iV- lO.STT/TV^). Error bars show one standard deviation (statistical error). 
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